Chapter 4 Diversity analysis
genome_counts_filt <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0)%>%
select(where(~ sum(.) > 0)) %>%
select(-AC85) %>%
rownames_to_column(., "genome")
genome_tree <- keep.tip(genome_tree, tip=genome_counts_filt$genome)
table <- genome_counts_filt %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character,as.numeric) %>%
rownames_to_column(., "sample")
master_table <- sample_metadata %>%
mutate(sample=Tube_code) %>%
mutate(Tube_code=str_remove_all(Tube_code, "_a")) %>%
filter(Tube_code %in% table$sample) %>%
mutate(treatment = sub("^\\d+_", "", time_point))%>%
left_join(., table, by=join_by("Tube_code" =="sample"))
sample_metadata <- master_table %>%
select(1:13)
genome_counts_filt <- master_table %>%
select(12,14:323) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character,as.numeric) %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0)%>%
dplyr::select(where(~ !all(. == 0)))%>%
rownames_to_column(., "genome")
genome_counts_filtering <- master_table %>%
select(12,14:323) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character,as.numeric) %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0)%>%
dplyr::select(where(~ !all(. == 0)))
genome_tree <- keep.tip(genome_tree, tip=genome_counts_filt$genome)
#match_data(genome_counts_filtering,genome_tree)
genome_metadata <- genome_metadata %>%
filter(genome %in% genome_counts_filt$genome)4.1 Calculate diversities
4.1.1 Alpha diversity
# Calculate Hill numbers
richness <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 0) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(richness = 1) %>%
rownames_to_column(var = "sample")
neutral <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 1) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(neutral = 1) %>%
rownames_to_column(var = "sample")
phylogenetic <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 1, tree = genome_tree) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(phylogenetic = 1) %>%
rownames_to_column(var = "sample")
genome_gifts1 <- genome_gifts[genome_counts_filt$genome,]
genome_gifts1 <- genome_gifts1[, colSums(genome_gifts1 != 0) > 0]
genome_counts_filt1 <- genome_counts_filt[genome_counts_filt$genome %in% rownames(genome_gifts1),]
genome_counts_filt1 <- genome_counts_filt1 %>%
remove_rownames() %>%
column_to_rownames(var = "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
rownames_to_column(., "genome")
# Merge all metrics
alpha_div <- richness %>%
full_join(neutral, by = join_by(sample == sample)) %>%
full_join(phylogenetic, by = join_by(sample == sample))4.1.2 Beta diversity
beta_q0n <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0)) %>%
hillpair(., q = 0)
beta_q1n <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0)) %>%
hillpair(., q = 1)
beta_q1p <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0)) %>%
hillpair(., q = 1, tree = genome_tree)
genome_gifts1 <- genome_gifts[genome_counts_filt$genome,]
genome_gifts1 <- genome_gifts1[, colSums(genome_gifts1 != 0) > 0]
genome_counts_filt1 <- genome_counts_filt[genome_counts_filt$genome %in% rownames(genome_gifts1),]4.2 Does captivity change the microbial community?
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == sample))%>%
filter(time_point %in% c("Acclimation", "Wild") ) 4.2.1 Alpha diversity
label_map <- c(
"Cold_wet" = "Cold population",
"Hot_dry" = "Warm-population",
"richness" = "Species Richness",
"neutral" = "Neutral Diversity",
"phylogenetic" = "Phylogenetic Diversity")
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
mutate(
metric = factor(metric, levels = c("richness", "neutral", "phylogenetic")),
time_point = factor(time_point, levels = c("Wild", "Acclimation"))) %>%
filter(time_point %in% c("Wild", "Acclimation")) %>%
ggplot(aes(y = value, x = time_point, color=Population, fill=Population)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#008080', "#d57d2c")) +
scale_fill_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#00808050', "#d57d2c50")) +
facet_grid(metric ~ Population, scales = "free_y", labeller = labeller(metric = label_map, type = label_map))+
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")
Richness
Modelq0GLMMNB <- glmer.nb(richness ~ Population*time_point+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Negative Binomial(9.4635) ( log )
Formula: richness ~ Population * time_point + (1 | individual)
Data: alpha_div_meta
AIC BIC logLik deviance df.resid
518.4 530.2 -253.2 506.4 47
Scaled residuals:
Min 1Q Median 3Q Max
-2.07994 -0.34567 0.00022 0.55317 1.22959
Random effects:
Groups Name Variance Std.Dev.
individual (Intercept) 0.1126 0.3355
Number of obs: 53, groups: individual, 27
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.7707 0.1208 31.218 < 2e-16 ***
PopulationHot_dry 0.6756 0.2006 3.368 0.000757 ***
time_pointWild 0.4312 0.1223 3.525 0.000424 ***
PopulationHot_dry:time_pointWild -0.5331 0.2048 -2.604 0.009221 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) PpltH_ tm_pnW
PpltnHt_dry -0.602
time_pntWld -0.557 0.335
PpltnHt_:_W 0.338 -0.526 -0.601$emmeans
time_point emmean SE df asymp.LCL asymp.UCL
Acclimation 4.11 0.100 Inf 3.91 4.31
Wild 4.27 0.099 Inf 4.08 4.47
Results are averaged over the levels of: Population
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
contrast estimate SE df z.ratio p.value
Acclimation - Wild -0.165 0.102 Inf -1.615 0.1063
Results are averaged over the levels of: Population
Results are given on the log (not the response) scale. $emmeans
Population emmean SE df asymp.LCL asymp.UCL
Cold_wet 3.99 0.101 Inf 3.79 4.18
Hot_dry 4.40 0.138 Inf 4.12 4.67
Results are averaged over the levels of: time_point
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
contrast estimate SE df z.ratio p.value
Cold_wet - Hot_dry -0.409 0.171 Inf -2.396 0.0166
Results are averaged over the levels of: time_point
Results are given on the log (not the response) scale. Neutral
Model_neutral <- nlme::lme(fixed = neutral ~ time_point*Population, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_neutral)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
402.7832 414.1342 -195.3916
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 5.51778 10.54924
Fixed effects: neutral ~ time_point * Population
Value Std.Error DF t-value p-value
(Intercept) 19.37729 2.883717 25 6.719554 0.0000
time_pointWild 17.27881 3.578682 24 4.828260 0.0001
PopulationHot_dry 24.71329 4.905493 25 5.037882 0.0000
time_pointWild:PopulationHot_dry -22.28917 6.126768 24 -3.637998 0.0013
Correlation:
(Intr) tm_pnW PpltH_
time_pointWild -0.642
PopulationHot_dry -0.588 0.377
time_pointWild:PopulationHot_dry 0.375 -0.584 -0.632
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.5974563 -0.5606149 0.1376820 0.5636829 1.6703831
Number of Observations: 53
Number of Groups: 27 $emmeans
time_point emmean SE df lower.CL upper.CL
Acclimation 31.7 2.45 25 26.7 36.8
Wild 37.9 2.43 24 32.9 42.9
Results are averaged over the levels of: Population
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Acclimation - Wild -6.13 3.06 24 -2.002 0.0567
Results are averaged over the levels of: Population
Degrees-of-freedom method: containment $emmeans
Population emmean SE df lower.CL upper.CL
Cold_wet 28.0 2.21 24 23.5 32.6
Hot_dry 41.6 3.09 24 35.2 48.0
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Cold_wet - Hot_dry -13.6 3.8 24 -3.568 0.0016
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment Phylogenetic
Model_phylo <- nlme::lme(fixed = phylogenetic ~ time_point*Population, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_phylo)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
193.7567 205.1077 -90.87837
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 0.7860836 1.191838
Fixed effects: phylogenetic ~ time_point * Population
Value Std.Error DF t-value p-value
(Intercept) 5.343318 0.3453897 25 15.470404 0.0000
time_pointWild -0.135493 0.4048212 24 -0.334698 0.7408
PopulationHot_dry 1.119702 0.5880336 25 1.904146 0.0685
time_pointWild:PopulationHot_dry -0.602016 0.6924897 24 -0.869350 0.3933
Correlation:
(Intr) tm_pnW PpltH_
time_pointWild -0.608
PopulationHot_dry -0.587 0.357
time_pointWild:PopulationHot_dry 0.355 -0.585 -0.596
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.2027374 -0.4577650 0.0610368 0.3504732 1.7675077
Number of Observations: 53
Number of Groups: 27 $emmeans
time_point emmean SE df lower.CL upper.CL
Acclimation 31.7 2.45 25 26.7 36.8
Wild 37.9 2.43 24 32.9 42.9
Results are averaged over the levels of: Population
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Acclimation - Wild -6.13 3.06 24 -2.002 0.0567
Results are averaged over the levels of: Population
Degrees-of-freedom method: containment $emmeans
Population emmean SE df lower.CL upper.CL
Cold_wet 28.0 2.21 24 23.5 32.6
Hot_dry 41.6 3.09 24 35.2 48.0
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Cold_wet - Hot_dry -13.6 3.8 24 -3.568 0.0016
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment 4.2.2 Beta diversity
samples_to_keep <- sample_metadata %>%
filter(time_point %in% c("Acclimation", "Wild")) %>%
select(Tube_code) %>%
pull()
subset_meta <- sample_metadata %>%
filter(time_point %in% c("Acclimation", "Wild"))Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$time_point) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.03445 0.034446 4.6041 999 0.036 *
Residuals 51 0.38156 0.007482
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(richness ~ time_point*Population,
data = subset_meta %>% arrange(match(Tube_code,labels(richness))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(richness))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| time_point | 1 | 0.5784094 | 0.03513500 | 2.236493 | 0.001 |
| Population | 1 | 2.8737192 | 0.17456169 | 11.111600 | 0.001 |
| time_point:Population | 1 | 0.3378127 | 0.02052015 | 1.306196 | 0.001 |
| Residual | 49 | 12.6725436 | 0.76978316 | NA | NA |
| Total | 52 | 16.4624849 | 1.00000000 | NA | NA |
subset_meta_arrange <- column_to_rownames(subset_meta, "Tube_code")
subset_meta_arrange<-subset_meta_arrange[labels(richness),]
pairwise <- pairwise.adonis(richness, subset_meta_arrange$Population, perm=999)
pairwise%>%
tt()| pairs | Df | SumsOfSqs | F.Model | R2 | p.value | p.adjusted | sig |
|---|---|---|---|---|---|---|---|
| Cold_wet vs Hot_dry | 1 | 2.87727 | 10.8015 | 0.1747774 | 0.001 | 0.001 | ** |
| pairs | Df | SumsOfSqs | F.Model | R2 | p.value | p.adjusted | sig |
|---|---|---|---|---|---|---|---|
| Wild vs Acclimation | 1 | 0.5784094 | 1.857135 | 0.035135 | 0.023 | 0.023 | . |
Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$time_point) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.01130 0.0113018 1.2646 999 0.259
Residuals 51 0.45578 0.0089369 adonis2(neutral ~ time_point*Population,
data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(neutral))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| time_point | 1 | 0.8535852 | 0.05306573 | 3.688514 | 0.001 |
| Population | 1 | 3.3711165 | 0.20957575 | 14.567274 | 0.001 |
| time_point:Population | 1 | 0.5212922 | 0.03240772 | 2.252609 | 0.001 |
| Residual | 49 | 11.3394385 | 0.70495080 | NA | NA |
| Total | 52 | 16.0854325 | 1.00000000 | NA | NA |
subset_meta_arrange <- column_to_rownames(subset_meta, "Tube_code")
subset_meta_arrange<-subset_meta_arrange[labels(neutral),]
pairwise <- pairwise.adonis(neutral, subset_meta_arrange$Population, perm=999)
pairwise%>%
tt()| pairs | Df | SumsOfSqs | F.Model | R2 | p.value | p.adjusted | sig |
|---|---|---|---|---|---|---|---|
| Cold_wet vs Hot_dry | 1 | 3.369473 | 13.51397 | 0.2094736 | 0.001 | 0.001 | ** |
| pairs | Df | SumsOfSqs | F.Model | R2 | p.value | p.adjusted | sig |
|---|---|---|---|---|---|---|---|
| Wild vs Acclimation | 1 | 0.8535852 | 2.858015 | 0.05306573 | 0.002 | 0.002 | * |
Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$time_point) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.00207 0.0020675 0.125 999 0.735
Residuals 51 0.84374 0.0165440 adonis2(phylo ~ time_point*Population,
data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(phylo))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| time_point | 1 | 0.2282721 | 0.07008479 | 4.486412 | 0.001 |
| Population | 1 | 0.3478468 | 0.10679697 | 6.836508 | 0.001 |
| time_point:Population | 1 | 0.1878081 | 0.05766139 | 3.691140 | 0.003 |
| Residual | 49 | 2.4931580 | 0.76545685 | NA | NA |
| Total | 52 | 3.2570850 | 1.00000000 | NA | NA |
subset_meta_arrange <- column_to_rownames(subset_meta, "Tube_code")
subset_meta_arrange<-subset_meta_arrange[labels(phylo),]
pairwise <- pairwise.adonis(phylo, subset_meta_arrange$Population, perm=999)
pairwise%>%
tt()| pairs | Df | SumsOfSqs | F.Model | R2 | p.value | p.adjusted | sig |
|---|---|---|---|---|---|---|---|
| Cold_wet vs Hot_dry | 1 | 0.3462509 | 6.066576 | 0.106307 | 0.002 | 0.002 | * |
| pairs | Df | SumsOfSqs | F.Model | R2 | p.value | p.adjusted | sig |
|---|---|---|---|---|---|---|---|
| Wild vs Acclimation | 1 | 0.2282721 | 3.84371 | 0.07008479 | 0.004 | 0.004 | * |
dbRDA
#Richness
cca_ord <- capscale(formula = richness ~ subset_meta$time_point* subset_meta$Population)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta$time_pointWild" = "Wild",
"subset_meta$PopulationHot_dry"="Warm population",
"subset_meta$time_pointWild:subset_meta$PopulationHot_dry"="Interaction Wam population")
beta_richness_rda <-CAP_df %>%
group_by(Population, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=Population,shape = time_point)) +
scale_color_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#008080', "#d57d2c")) +
scale_fill_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#00808050', "#d57d2c50")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()
#Neutral
cca_ord <- capscale(formula = neutral ~ subset_meta$time_point* subset_meta$Population)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta$time_pointWild" = "Wild",
"subset_meta$PopulationHot_dry"="Warm population",
"subset_meta$time_pointWild:subset_meta$PopulationHot_dry"="Interaction Wam population")
beta_neutral_rda <- CAP_df %>%
group_by(Population, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=Population,shape = time_point)) +
scale_color_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#008080', "#d57d2c")) +
scale_fill_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#00808050', "#d57d2c50"))+
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()
#Phylogenetic
cca_ord <- capscale(formula = phylo ~ subset_meta$time_point* subset_meta$Population)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta$time_pointWild" = "Wild",
"subset_meta$PopulationHot_dry"="Warm population",
"subset_meta$time_pointWild:subset_meta$PopulationHot_dry"="Interaction Wam population")
beta_phylogenetic_rda<- CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=Population,shape = time_point)) +
scale_color_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#008080', "#d57d2c")) +
scale_fill_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#00808050', "#d57d2c50")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()ggarrange(beta_richness_rda, beta_neutral_rda, beta_phylogenetic_rda, ncol=3, nrow=1, common.legend = TRUE, legend="right")
4.2.3 Cold population
4.2.3.1 Structural zeros
wild_samples <- sample_metadata %>%
filter(Population == "Cold_wet" & time_point == "Wild") %>%
dplyr::select(sample) %>%
pull()
acclimation_samples <- sample_metadata %>%
filter(Population == "Cold_wet" & time_point == "Acclimation") %>%
dplyr::select(sample) %>% pull()
structural_zeros <- genome_counts_filt %>%
select(one_of(c("genome",subset_meta$sample))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
rowwise() %>% #compute for each row (genome)
mutate(all_zeros_wild = all(c_across(all_of(wild_samples)) == 0)) %>%
mutate(all_zeros_acclimation= all(c_across(all_of(acclimation_samples)) == 0)) %>%
mutate(average_wild = mean(c_across(all_of(wild_samples)), na.rm = TRUE)) %>%
mutate(average_acclimation = mean(c_across(all_of(acclimation_samples)), na.rm = TRUE)) %>%
filter(all_zeros_wild == TRUE || all_zeros_acclimation==TRUE) %>%
mutate(present = case_when(
all_zeros_wild & !all_zeros_acclimation ~ "acclimation",
!all_zeros_wild & all_zeros_acclimation ~ "wild",
!all_zeros_wild & !all_zeros_acclimation ~ "None",
TRUE ~ NA_character_
)) %>%
mutate(average = ifelse(present == "acclimation", average_wild, average_acclimation)) %>%
dplyr::select(genome, present, average) %>%
left_join(genome_metadata, by=join_by(genome==genome)) %>%
arrange(present,-average)Wild
structural_zeros %>%
filter(present=="wild") %>%
count(phylum, name = "sample_count") %>%
arrange(desc(sample_count)) # A tibble: 5 × 2
# Rowwise:
phylum sample_count
<chr> <int>
1 p__Bacillota_A 7
2 p__Bacillota 1
3 p__Bacillota_B 1
4 p__Bacteroidota 1
5 p__Spirochaetota 1Acclimation
# A tibble: 14 × 7
# Rowwise:
genome phylum class order family genus species
<chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 AH1_2nd_12:bin_000001 p__Pseudomonadota c__Gammaproteobacteria o__Pseudomonadales f__Moraxellaceae g__Acinetobacter s__Acinetobacter johnsonii
2 AH1_2nd_15:bin_000001 p__Pseudomonadota c__Alphaproteobacteria o__Rhizobiales f__Rhizobiaceae g__Agrobacterium s__Agrobacterium tumefaciens_H
3 AH1_2nd_17:bin_000010 p__Bacteroidota c__Bacteroidia o__Bacteroidales f__Rikenellaceae g__Mucinivorans s__
4 AH1_2nd_17:bin_000038 p__Bacteroidota c__Bacteroidia o__Bacteroidales f__Bacteroidaceae g__Bacteroides_G s__
5 AH1_2nd_20:bin_000014 p__Pseudomonadota c__Gammaproteobacteria o__Enterobacterales f__Enterobacteriaceae g__Citrobacter s__Citrobacter portucalensis
6 AH1_2nd_20:bin_000061 p__Bacillota c__Bacilli o__Lactobacillales f__Enterococcaceae g__Enterococcus s__
7 AH1_2nd_20:bin_000073 p__Bacillota c__Bacilli o__Lactobacillales f__Enterococcaceae g__Enterococcus s__Enterococcus sp002174455
8 LI1_2nd_10:bin_000001 p__Bacillota c__Bacilli o__Lactobacillales f__Enterococcaceae g__Enterococcus_A s__Enterococcus_A raffinosus
9 LI1_2nd_10:bin_000017 p__Chlamydiota c__Chlamydiia o__Chlamydiales f__Chlamydiaceae g__ s__
10 LI1_2nd_2:bin_000039 p__Bacillota_A c__Clostridia o__TANB77 f__CAG-508 g__ s__
11 LI1_2nd_7:bin_000045 p__Bacillota_A c__Clostridia o__Oscillospirales f__Acutalibacteraceae g__ s__
12 LI1_2nd_7:bin_000074 p__Pseudomonadota c__Gammaproteobacteria o__Enterobacterales f__Enterobacteriaceae g__Escherichia s__Escherichia coli
13 LI1_2nd_8:bin_000030 p__Actinomycetota c__Actinomycetia o__Mycobacteriales f__Mycobacteriaceae g__Corynebacterium s__
14 LI1_2nd_8:bin_000079 p__Pseudomonadota c__Gammaproteobacteria o__Enterobacterales f__Enterobacteriaceae g__Citrobacter_A s__Citrobacter_A amalonaticus Community elements differences in structural zeros
element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
element_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)# A tibble: 0 × 3
# ℹ 3 variables: trait <chr>, p_value <dbl>, p_adjust <dbl>4.2.3.2 Differential abundance analysis: composition
phylo_samples <- sample_metadata %>%
filter(Population == "Cold_wet" & time_point %in% c("Wild", "Acclimation") )%>%
column_to_rownames("sample") %>%
sample_data()
phylo_genome <- genome_counts_filt %>%
filter(!genome %in% structural_zeros$genome) %>%
select(one_of(c("genome",rownames(phylo_samples)))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
column_to_rownames("genome") %>%
otu_table(., taxa_are_rows = TRUE)
phylo_taxonomy <- genome_metadata %>%
filter(genome %in% rownames(phylo_genome)) %>%
column_to_rownames("genome") %>%
dplyr::select(domain,phylum,class,order,family,genus,species) %>%
as.matrix() %>%
tax_table()
physeq_genome_filtered <- phyloseq(phylo_genome, phylo_taxonomy, phylo_samples)ancom_rand_output = ancombc2(data = physeq_genome_filtered,
assay_name = "counts",
tax_level = NULL,
fix_formula = "time_point",
rand_formula = "(1|individual)",
p_adj_method = "holm",
pseudo_sens = TRUE,
prv_cut =0,
lib_cut = 0,
s0_perc = 0.05,
group = NULL,
struc_zero = FALSE,
neg_lb = FALSE,
alpha = 0.05,
n_cl = 2,
verbose = TRUE,
global = FALSE,
pairwise = FALSE,
dunnet = FALSE,
trend = FALSE,
iter_control = list(tol = 1e-5, max_iter = 20, verbose = FALSE),
em_control = list(tol = 1e-5, max_iter = 100),
lme_control = lme4::lmerControl(),
mdfdr_control = list(fwer_ctrl_method = "holm", B = 100),
trend_control = NULL)taxonomy <- data.frame(physeq_genome_filtered@tax_table) %>%
rownames_to_column(., "taxon") %>%
mutate_at(vars(order, phylum, family, genus, species), ~ str_replace(., "[dpcofgs]__", ""))
ancombc_rand_table_mag <- ancom_rand_output$res %>%
dplyr::select(taxon, lfc_time_point1_Acclimation, p_time_point1_Acclimation) %>%
filter(p_time_point1_Acclimation < 0.05) %>%
dplyr::arrange(p_time_point1_Acclimation) %>%
merge(., taxonomy, by="taxon") %>%
mutate_at(vars(phylum, species), ~ str_replace(., "[dpcofgs]__", ""))%>%
dplyr::arrange(lfc_time_point1_Acclimation)
colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
mutate_at(vars(phylum), ~ str_replace(., "[dpcofgs]__", "")) %>%
right_join(taxonomy, by=join_by(phylum == phylum)) %>%
dplyr::select(phylum, colors) %>%
mutate(colors = str_c(colors, "80")) %>% #add 80% alpha
unique() %>%
dplyr::arrange(phylum)
tax_table <- as.data.frame(unique(ancombc_rand_table_mag$phylum))
colnames(tax_table)[1] <- "phylum"
tax_color <- merge(tax_table, colors_alphabetic, by="phylum")%>%
dplyr::arrange(phylum) %>%
dplyr::select(colors) %>%
pull()ancombc_rand_table_mag%>%
mutate(genome=factor(taxon,levels=ancombc_rand_table_mag$taxon)) %>%
ggplot(., aes(x=lfc_time_point1_Acclimation, y=forcats::fct_reorder(genome,lfc_time_point1_Acclimation), fill=phylum)) + #forcats::fct_rev()
geom_col() +
scale_fill_manual(values=tax_color) +
geom_hline(yintercept=0) +
# coord_flip()+
theme(panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 14, face = "bold"),
legend.text = element_text(size = 12),
legend.title = element_text(size = 14, face = "bold"),
legend.position = "right", legend.box = "vertical")+
xlab("log2FoldChange") +
ylab("Species")+
guides(fill=guide_legend(title="Phylum"))
Phyla of the significant MAGs
ancombc_rand_table_mag%>%
filter(lfc_time_point1_Acclimation<0) %>%
count(phylum, name = "sample_count") %>%
arrange(desc(sample_count)) phylum sample_count
1 Bacillota_A 16
2 Bacteroidota 9
3 Bacillota_C 1 phylum genus species
1 Bacillota_A MGBC140009
2 Bacillota_A
3 Bacillota_A
4 Bacillota_A Negativibacillus
5 Bacteroidota Parabacteroides
6 Bacteroidota Bacteroides
7 Bacillota_A Hungatella
8 Bacteroidota Bacteroides
9 Bacteroidota Parabacteroides
10 Bacteroidota Parabacteroides
11 Bacteroidota Bacteroides
12 Bacillota_A
13 Bacteroidota Bacteroides
14 Bacillota_A
15 Bacillota_A Copromonas
16 Bacillota_A Enterocloster
17 Bacillota_A Velocimicrobium
18 Bacillota_A CAG-95
19 Bacillota_C
20 Bacteroidota Bacteroides
21 Bacillota_A
22 Bacillota_A Copromonas
23 Bacteroidota Bacteroides Bacteroides thetaiotaomicron
24 Bacillota_A Pseudoflavonifractor
25 Bacillota_A Enterocloster
26 Bacillota_A JALFVM01 ancombc_rand_table_mag%>%
filter(lfc_time_point1_Acclimation<0) %>%
count(genus, name = "sample_count") %>%
arrange(desc(sample_count)) genus sample_count
1 6
2 Bacteroides 6
3 Parabacteroides 3
4 Copromonas 2
5 Enterocloster 2
6 CAG-95 1
7 Hungatella 1
8 JALFVM01 1
9 MGBC140009 1
10 Negativibacillus 1
11 Pseudoflavonifractor 1
12 Velocimicrobium 14.2.3.3 Differences in the functional capacity
GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts_filt$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>%
select_if(~ !is.numeric(.) || sum(.) != 0)
GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)
GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)
sample_sub <- sample_metadata %>%
filter(Population == "Cold_wet" & treatment %in% c("Acclimation", "Wild"))
genome_counts_row <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>%
select(one_of(c("genome",sample_sub$sample))) %>%
column_to_rownames(., "genome")
GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_row,GIFT_db)4.2.3.3.1 Community elements differences:
MCI
GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
treatment MCI sd
<chr> <dbl> <dbl>
1 Acclimation 0.335 0.0324
2 Wild 0.350 0.0237MCI <- GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.94347, p-value = 0.07144
Wilcoxon rank sum exact test
data: value by treatment
W = 114, p-value = 0.2066
alternative hypothesis: true location shift is not equal to 0element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
significant_elements <- element_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))
element_gift_t <- element_gift %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "trait")
element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "sample")%>%
left_join(., sample_metadata[c(12,13)], by = join_by(sample == sample))
difference_table <- element_gift_filt %>%
select(-sample) %>%
group_by(treatment) %>%
summarise(across(everything(), mean)) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>%
arrange(Function) %>%
mutate(Difference=Acclimation-Wild)%>%
mutate(group_color = ifelse(Difference <0, "Wild","Acclimation")) difference_table %>%
ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) +
geom_col() +
# geom_point(size=4) +
scale_fill_manual(values=c('#008080', '#00808050')) +
geom_hline(yintercept=0) +
coord_flip()+
theme(axis.text = element_text(size = 10),
axis.title = element_text(size = 12),
legend.position = "right",
panel.background = element_blank(),
panel.grid.major = element_line(size = 0.15, linetype = 'solid',
colour = "grey"))+
xlab("Microbial Functional Capacity") +
ylab("Mean difference")+
labs(fill="Treatment")
4.2.3.3.2 Community functions differences
MCI
GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
treatment MCI sd
<chr> <dbl> <dbl>
1 Acclimation 0.330 0.0320
2 Wild 0.346 0.0196MCI <- GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.94507, p-value = 0.07998
Wilcoxon rank sum exact test
data: value by treatment
W = 100, p-value = 0.08313
alternative hypothesis: true location shift is not equal to 0function_gift <- GIFTs_functions_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
merge(., sample_metadata[c(12,13)], by="sample")unique_funct_db<- GIFT_db[c(3,4,5)] %>%
distinct(Code_function, .keep_all = TRUE)
significant_functional <- function_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
left_join(.,unique_funct_db[c(1,3)],by = join_by(trait == Code_function))| Code_function | Acclimation | Wild | Function | Difference | treatment |
|---|---|---|---|---|---|
| B03 | 0.1097356 | 0.1256091 | Amino acid derivative biosynthesis | -0.0158735 | Wild |
| D03 | 0.2921075 | 0.3442827 | Sugar degradation | -0.0521752 | Wild |
4.2.4 Warm population
4.2.4.1 Structural zeros
wild_samples_warm <- sample_metadata %>%
filter(Population == "Hot_dry" & time_point == "Wild") %>%
dplyr::select(sample) %>%
pull()
acclimation_samples_warm <- sample_metadata %>%
filter(Population == "Hot_dry" & time_point == "Acclimation") %>%
dplyr::select(sample) %>% pull()
structural_zeros <- genome_counts_filt %>%
select(one_of(c("genome",subset_meta$sample))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
rowwise() %>% #compute for each row (genome)
mutate(all_zeros_wild_warm = all(c_across(all_of(wild_samples_warm)) == 0)) %>%
mutate(all_zeros_acclimation_warm= all(c_across(all_of(acclimation_samples_warm)) == 0)) %>%
mutate(average_wild_warm = mean(c_across(all_of(wild_samples_warm)), na.rm = TRUE)) %>%
mutate(average_acclimation_warm = mean(c_across(all_of(acclimation_samples_warm)), na.rm = TRUE)) %>%
filter(all_zeros_wild_warm == TRUE || all_zeros_acclimation_warm==TRUE) %>%
mutate(present = case_when(
all_zeros_wild_warm & !all_zeros_acclimation_warm ~ "acclimation",
!all_zeros_wild_warm & all_zeros_acclimation_warm ~ "wild",
!all_zeros_wild_warm & !all_zeros_acclimation_warm ~ "None",
TRUE ~ NA_character_
)) %>%
mutate(average = ifelse(present == "acclimation", average_wild_warm, average_acclimation_warm)) %>%
dplyr::select(genome, present, average) %>%
left_join(genome_metadata, by=join_by(genome==genome)) %>%
arrange(present,-average)Wild
structural_zeros %>%
filter(present=="wild") %>%
count(phylum, name = "sample_count") %>%
arrange(desc(sample_count)) # A tibble: 7 × 2
# Rowwise:
phylum sample_count
<chr> <int>
1 p__Bacillota_A 11
2 p__Bacillota 4
3 p__Pseudomonadota 4
4 p__Bacillota_C 2
5 p__Bacteroidota 1
6 p__Desulfobacterota 1
7 p__Fusobacteriota 1Acclimation
# A tibble: 32 × 7
# Rowwise:
genome phylum class order family genus species
<chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 AH1_2nd_10:bin_000010 p__Cyanobacteriota c__Vampirovibrionia o__Gastranaerophilales f__Gastranaerophilaceae g__CALURL01 s__
2 AH1_2nd_12:bin_000001 p__Pseudomonadota c__Gammaproteobacteria o__Pseudomonadales f__Moraxellaceae g__Acinetobacter s__Acinetobacter johnsonii
3 AH1_2nd_14:bin_000015 p__Bacillota_A c__Clostridia o__Christensenellales f__UBA3700 g__ s__
4 AH1_2nd_14:bin_000063 p__Bacillota_A c__Clostridia o__Lachnospirales f__Lachnospiraceae g__Blautia_A s__
5 AH1_2nd_15:bin_000001 p__Pseudomonadota c__Alphaproteobacteria o__Rhizobiales f__Rhizobiaceae g__Agrobacterium s__Agrobacterium tumefaciens_H
6 AH1_2nd_16:bin_000033 p__Bacteroidota c__Bacteroidia o__Bacteroidales f__ g__ s__
7 AH1_2nd_16:bin_000096 p__Bacillota_A c__Clostridia o__Oscillospirales f__Ruminococcaceae g__Ruthenibacterium s__
8 AH1_2nd_17:bin_000015 p__Bacillota c__Bacilli o__Erysipelotrichales f__Erysipelotrichaceae g__ s__
9 AH1_2nd_17:bin_000023 p__Bacillota c__Bacilli o__Erysipelotrichales f__Coprobacillaceae g__Coprobacillus s__
10 AH1_2nd_18:bin_000039 p__Bacteroidota c__Bacteroidia o__Bacteroidales f__Rikenellaceae g__Alistipes s__
# ℹ 22 more rowsCommunity elements differences in structural zeros
element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
element_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)# A tibble: 0 × 3
# ℹ 3 variables: trait <chr>, p_value <dbl>, p_adjust <dbl>4.2.5 Differential abundance analysis: composition
phylo_samples_warm <- sample_metadata %>%
filter(Population == "Hot_dry" & time_point %in% c("Wild", "Acclimation") )%>%
column_to_rownames("sample") %>%
sample_data()
phylo_genome_warm <- genome_counts_filt %>%
filter(!genome %in% structural_zeros$genome) %>%
select(one_of(c("genome",rownames(phylo_samples_warm)))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
column_to_rownames("genome") %>%
otu_table(., taxa_are_rows = TRUE)
phylo_taxonomy_warm <- genome_metadata %>%
filter(genome %in% rownames(phylo_genome_warm)) %>%
column_to_rownames("genome") %>%
dplyr::select(domain,phylum,class,order,family,genus,species) %>%
as.matrix() %>%
tax_table()
physeq_genome_filtered_warm <- phyloseq(phylo_genome_warm, phylo_taxonomy_warm, phylo_samples_warm)ancom_rand_output_warm = ancombc2(data = physeq_genome_filtered_warm,
assay_name = "counts",
tax_level = NULL,
fix_formula = "time_point",
rand_formula = "(1|individual)",
p_adj_method = "holm",
pseudo_sens = TRUE,
prv_cut =0,
lib_cut = 0,
s0_perc = 0.05,
group = NULL,
struc_zero = FALSE,
neg_lb = FALSE,
alpha = 0.05,
n_cl = 2,
verbose = TRUE,
global = FALSE,
pairwise = FALSE,
dunnet = FALSE,
trend = FALSE,
iter_control = list(tol = 1e-5, max_iter = 20, verbose = FALSE),
em_control = list(tol = 1e-5, max_iter = 100),
lme_control = lme4::lmerControl(),
mdfdr_control = list(fwer_ctrl_method = "holm", B = 100),
trend_control = NULL)taxonomy <- data.frame(physeq_genome_filtered_warm@tax_table) %>%
rownames_to_column(., "taxon") %>%
mutate_at(vars(order, phylum, family, genus, species), ~ str_replace(., "[dpcofgs]__", ""))
ancombc_rand_table_mag <- ancom_rand_output_warm$res %>%
dplyr::select(taxon, lfc_time_point1_Acclimation, p_time_point1_Acclimation) %>%
filter(p_time_point1_Acclimation < 0.05) %>%
dplyr::arrange(p_time_point1_Acclimation) %>%
merge(., taxonomy, by="taxon") %>%
mutate_at(vars(phylum, species), ~ str_replace(., "[dpcofgs]__", ""))%>%
dplyr::arrange(lfc_time_point1_Acclimation)
colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
mutate_at(vars(phylum), ~ str_replace(., "[dpcofgs]__", "")) %>%
right_join(taxonomy, by=join_by(phylum == phylum)) %>%
dplyr::select(phylum, colors) %>%
mutate(colors = str_c(colors, "80")) %>% #add 80% alpha
unique() %>%
dplyr::arrange(phylum)
tax_table <- as.data.frame(unique(ancombc_rand_table_mag$phylum))
colnames(tax_table)[1] <- "phylum"
tax_color <- merge(tax_table, colors_alphabetic, by="phylum")%>%
dplyr::arrange(phylum) %>%
dplyr::select(colors) %>%
pull()ancombc_rand_table_mag%>%
mutate(genome=factor(taxon,levels=ancombc_rand_table_mag$taxon)) %>%
ggplot(., aes(x=lfc_time_point1_Acclimation, y=forcats::fct_reorder(genome,lfc_time_point1_Acclimation), fill=phylum)) + #forcats::fct_rev()
geom_col() +
scale_fill_manual(values=tax_color) +
geom_hline(yintercept=0) +
# coord_flip()+
theme(panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 14, face = "bold"),
legend.text = element_text(size = 12),
legend.title = element_text(size = 14, face = "bold"),
legend.position = "right", legend.box = "vertical")+
xlab("log2FoldChange") +
ylab("Species")+
guides(fill=guide_legend(title="Phylum"))Phyla of the significant MAGs
ancombc_rand_table_mag%>%
filter(lfc_time_point1_Acclimation<0) %>%
count(phylum, name = "sample_count") %>%
arrange(desc(sample_count))
ancombc_rand_table_mag%>%
filter(lfc_time_point1_Acclimation<0) %>%
select(phylum, genus, species)
ancombc_rand_table_mag%>%
filter(lfc_time_point1_Acclimation<0) %>%
count(genus, name = "sample_count") %>%
arrange(desc(sample_count))4.2.6 Differences in the functional capacity
GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts_filt$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>%
select_if(~ !is.numeric(.) || sum(.) != 0)
GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)
GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)
sample_sub <- sample_metadata %>%
filter(Population == "Hot_dry" & treatment %in% c("Acclimation", "Wild"))
genome_counts_row <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>%
select(one_of(c("genome",sample_sub$sample))) %>%
column_to_rownames(., "genome")
GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_row,GIFT_db)4.2.6.1 Community elements differences:
MCI
GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
treatment MCI sd
<chr> <dbl> <dbl>
1 Acclimation 0.355 0.0228
2 Wild 0.317 0.0335MCI <- GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.94062, p-value = 0.297
Wilcoxon rank sum exact test
data: value by treatment
W = 68, p-value = 0.01419
alternative hypothesis: true location shift is not equal to 0element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
significant_elements <- element_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))
element_gift_t <- element_gift %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "trait")
element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "sample")%>%
left_join(., sample_metadata[c(12,13)], by = join_by(sample == sample))
difference_table <- element_gift_filt %>%
select(-sample) %>%
group_by(treatment) %>%
summarise(across(everything(), mean)) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>%
arrange(Function) %>%
mutate(Difference=Acclimation-Wild)%>%
mutate(group_color = ifelse(Difference <0, "Wild","Acclimation")) difference_table %>%
ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) +
geom_col() +
# geom_point(size=4) +
scale_fill_manual(values=c("#d57d2c50","#d57d2c")) +
geom_hline(yintercept=0) +
coord_flip()+
theme(axis.text = element_text(size = 10),
axis.title = element_text(size = 12),
legend.position = "right",
panel.background = element_blank(),
panel.grid.major = element_line(size = 0.15, linetype = 'solid',
colour = "grey"))+
xlab("Microbial Functional Capacity") +
ylab("Mean difference")+
labs(fill="Treatment")
4.2.6.2 Community functions differences
MCI
GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
treatment MCI sd
<chr> <dbl> <dbl>
1 Acclimation 0.348 0.0158
2 Wild 0.327 0.0245MCI <- GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.95485, p-value = 0.5061
Wilcoxon rank sum exact test
data: value by treatment
W = 61, p-value = 0.07701
alternative hypothesis: true location shift is not equal to 0function_gift <- GIFTs_functions_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
merge(., sample_metadata[c(12,13)], by="sample")unique_funct_db<- GIFT_db[c(3,4,5)] %>%
distinct(Code_function, .keep_all = TRUE)
significant_functional <- function_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
left_join(.,unique_funct_db[c(1,3)],by = join_by(trait == Code_function))| Code_function | Acclimation | Wild | Function | Difference | treatment |
|---|---|---|---|---|---|
| S02 | 0.1699218 | 0.2472264 | Appendages | -0.0773046 | Wild |
4.3 Does the antibiotic administration change the microbial community?
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == sample))%>%
filter(time_point %in% c("Acclimation", "Antibiotics") ) 4.3.1 Alpha diversity
label_map <- c(
"Cold_wet" = "Cold population",
"Hot_dry" = "Warm-population",
"richness" = "Species Richness",
"neutral" = "Neutral Diversity",
"phylogenetic" = "Phylogenetic Diversity")
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic"))) %>%
filter(time_point %in% c("Acclimation", "Antibiotics") ) %>%
ggplot(aes(y = value, x = time_point, color=Population, fill=Population)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#008080', "#d57d2c")) +
scale_fill_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#00808050', "#d57d2c50")) +
facet_grid(metric ~ Population, scales = "free_y", labeller = labeller(metric = label_map, type = label_map))+
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")
Richness
Modelq0GLMMNB <- glmer.nb(richness ~ time_point*Population+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Negative Binomial(5.1645) ( log )
Formula: richness ~ time_point * Population + (1 | individual)
Data: alpha_div_meta
AIC BIC logLik deviance df.resid
434.8 446.1 -211.4 422.8 43
Scaled residuals:
Min 1Q Median 3Q Max
-1.61979 -0.56608 0.06074 0.51882 2.39923
Random effects:
Groups Name Variance Std.Dev.
individual (Intercept) 0.2193 0.4683
Number of obs: 49, groups: individual, 27
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.7459 0.1623 23.073 < 2e-16 ***
time_pointAntibiotics -1.1523 0.1862 -6.190 6.02e-10 ***
PopulationHot_dry 0.7055 0.2719 2.595 0.00946 **
time_pointAntibiotics:PopulationHot_dry -0.3586 0.3037 -1.181 0.23769
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) tm_pnA PpltH_
tm_pntAntbt -0.466
PpltnHt_dry -0.598 0.278
tm_pntA:PH_ 0.292 -0.611 -0.461$emmeans
time_point emmean SE df asymp.LCL asymp.UCL
Acclimation 4.10 0.136 Inf 3.83 4.36
Antibiotics 2.77 0.151 Inf 2.47 3.06
Results are averaged over the levels of: Population
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
contrast estimate SE df z.ratio p.value
Acclimation - Antibiotics 1.33 0.152 Inf 8.751 <.0001
Results are averaged over the levels of: Population
Results are given on the log (not the response) scale. $emmeans
Population emmean SE df asymp.LCL asymp.UCL
Cold_wet 3.17 0.145 Inf 2.89 3.45
Hot_dry 3.70 0.196 Inf 3.31 4.08
Results are averaged over the levels of: time_point
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
contrast estimate SE df z.ratio p.value
Cold_wet - Hot_dry -0.526 0.243 Inf -2.168 0.0302
Results are averaged over the levels of: time_point
Results are given on the log (not the response) scale. Neutral
Model_neutral <- nlme::lme(fixed = neutral ~ time_point*Population, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_neutral)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
342.5018 353.3418 -165.2509
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 3.282619 7.923552
Fixed effects: neutral ~ time_point * Population
Value Std.Error DF t-value p-value
(Intercept) 19.11713 2.078908 25 9.195758 0e+00
time_pointAntibiotics -11.46317 2.832826 20 -4.046551 6e-04
PopulationHot_dry 24.97345 3.534826 25 7.064972 0e+00
time_pointAntibiotics:PopulationHot_dry -20.91701 4.793363 20 -4.363745 3e-04
Correlation:
(Intr) tm_pnA PpltH_
time_pointAntibiotics -0.633
PopulationHot_dry -0.588 0.373
time_pointAntibiotics:PopulationHot_dry 0.374 -0.591 -0.632
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.8258256 -0.5685215 -0.1867832 0.5823712 2.0651455
Number of Observations: 49
Number of Groups: 27 $emmeans
time_point emmean SE df lower.CL upper.CL
Acclimation 31.60 1.77 25 27.96 35.2
Antibiotics 9.68 1.87 20 5.77 13.6
Results are averaged over the levels of: Population
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Acclimation - Antibiotics 21.9 2.4 20 9.147 <.0001
Results are averaged over the levels of: Population
Degrees-of-freedom method: containment $emmeans
Population emmean SE df lower.CL upper.CL
Cold_wet 13.4 1.61 20 10.0 16.7
Hot_dry 27.9 2.22 20 23.3 32.5
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Cold_wet - Hot_dry -14.5 2.74 20 -5.288 <.0001
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment Phylogenetic
Model_phylo <- nlme::lme(fixed = phylogenetic ~ time_point*Population, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_phylo)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
203.6546 214.4946 -95.8273
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 0.7173285 1.688356
Fixed effects: phylogenetic ~ time_point * Population
Value Std.Error DF t-value p-value
(Intercept) 5.365728 0.4446271 25 12.067930 0.0000
time_pointAntibiotics -0.761455 0.6038652 20 -1.260969 0.2218
PopulationHot_dry 1.097291 0.7560383 25 1.451370 0.1591
time_pointAntibiotics:PopulationHot_dry -0.702255 1.0216421 20 -0.687379 0.4997
Correlation:
(Intr) tm_pnA PpltH_
time_pointAntibiotics -0.631
PopulationHot_dry -0.588 0.371
time_pointAntibiotics:PopulationHot_dry 0.373 -0.591 -0.629
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.85993684 -0.51436620 -0.05130501 0.48944142 1.80209573
Number of Observations: 49
Number of Groups: 27 $emmeans
time_point emmean SE df lower.CL upper.CL
Acclimation 31.60 1.77 25 27.96 35.2
Antibiotics 9.68 1.87 20 5.77 13.6
Results are averaged over the levels of: Population
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Acclimation - Antibiotics 21.9 2.4 20 9.147 <.0001
Results are averaged over the levels of: Population
Degrees-of-freedom method: containment $emmeans
Population emmean SE df lower.CL upper.CL
Cold_wet 13.4 1.61 20 10.0 16.7
Hot_dry 27.9 2.22 20 23.3 32.5
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Cold_wet - Hot_dry -14.5 2.74 20 -5.288 <.0001
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment 4.3.2 Beta diversity
samples_to_keep <- sample_metadata %>%
filter(time_point %in% c("Acclimation", "Antibiotics")) %>%
select(Tube_code) %>%
pull()
subset_meta <- sample_metadata %>%
filter(time_point %in% c("Acclimation", "Antibiotics"))Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$time_point) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.03530 0.035300 10.073 999 0.004 **
Residuals 47 0.16471 0.003505
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(richness ~ time_point*Population,
data = subset_meta %>% arrange(match(Tube_code,labels(richness))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(richness))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| time_point | 1 | 1.9348458 | 0.1002347 | 6.065418 | 0.001 |
| Population | 1 | 2.1356198 | 0.1106358 | 6.694811 | 0.001 |
| time_point:Population | 1 | 0.8778553 | 0.0454773 | 2.751930 | 0.004 |
| Residual | 45 | 14.3548318 | 0.7436522 | NA | NA |
| Total | 48 | 19.3031527 | 1.0000000 | NA | NA |
subset_meta_arrange <- column_to_rownames(subset_meta, "Tube_code")
subset_meta_arrange<-subset_meta_arrange[labels(richness),]
pairwise <- pairwise.adonis(richness, subset_meta_arrange$Population, perm=999)
pairwise%>%
tt()| pairs | Df | SumsOfSqs | F.Model | R2 | p.value | p.adjusted | sig |
|---|---|---|---|---|---|---|---|
| Cold_wet vs Hot_dry | 1 | 2.135115 | 5.845187 | 0.1106096 | 0.001 | 0.001 | ** |
| pairs | Df | SumsOfSqs | F.Model | R2 | p.value | p.adjusted | sig |
|---|---|---|---|---|---|---|---|
| Acclimation vs Antibiotics | 1 | 1.934846 | 5.235844 | 0.1002347 | 0.001 | 0.001 | ** |
Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$time_point) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.051657 0.051657 9.9224 999 0.003 **
Residuals 47 0.244689 0.005206
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(neutral ~ time_point*Population,
data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(neutral))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| time_point | 1 | 2.0429592 | 0.11061666 | 7.254623 | 0.001 |
| Population | 1 | 2.8764490 | 0.15574623 | 10.214376 | 0.001 |
| time_point:Population | 1 | 0.8770558 | 0.04748846 | 3.114457 | 0.002 |
| Residual | 45 | 12.6723552 | 0.68614864 | NA | NA |
| Total | 48 | 18.4688192 | 1.00000000 | NA | NA |
subset_meta_arrange <- column_to_rownames(subset_meta, "Tube_code")
subset_meta_arrange<-subset_meta_arrange[labels(neutral),]
pairwise <- pairwise.adonis(neutral, subset_meta_arrange$Population, perm=999)
pairwise%>%
tt()| pairs | Df | SumsOfSqs | F.Model | R2 | p.value | p.adjusted | sig |
|---|---|---|---|---|---|---|---|
| Cold_wet vs Hot_dry | 1 | 2.876596 | 8.670991 | 0.1557542 | 0.001 | 0.001 | ** |
| pairs | Df | SumsOfSqs | F.Model | R2 | p.value | p.adjusted | sig |
|---|---|---|---|---|---|---|---|
| Acclimation vs Antibiotics | 1 | 2.042959 | 5.845605 | 0.1106167 | 0.001 | 0.001 | ** |
Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$time_point) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.67439 0.67439 55.932 999 0.001 ***
Residuals 47 0.56670 0.01206
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(phylo ~ time_point*Population,
data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(phylo))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| time_point | 1 | 1.8783556 | 0.23032747 | 16.358708 | 0.001 |
| Population | 1 | 0.7551942 | 0.09260332 | 6.577030 | 0.001 |
| time_point:Population | 1 | 0.3545686 | 0.04347786 | 3.087959 | 0.022 |
| Residual | 45 | 5.1670341 | 0.63359135 | NA | NA |
| Total | 48 | 8.1551525 | 1.00000000 | NA | NA |
subset_meta_arrange <- column_to_rownames(subset_meta, "Tube_code")
subset_meta_arrange<-subset_meta_arrange[labels(phylo),]
pairwise <- pairwise.adonis(phylo, subset_meta_arrange$Population, perm=999)
pairwise%>%
tt()| pairs | Df | SumsOfSqs | F.Model | R2 | p.value | p.adjusted | sig |
|---|---|---|---|---|---|---|---|
| Cold_wet vs Hot_dry | 1 | 0.753589 | 4.785297 | 0.09240649 | 0.001 | 0.001 | ** |
| pairs | Df | SumsOfSqs | F.Model | R2 | p.value | p.adjusted | sig |
|---|---|---|---|---|---|---|---|
| Acclimation vs Antibiotics | 1 | 1.878356 | 14.06493 | 0.2303275 | 0.001 | 0.001 | ** |
dbRDA
#Richness
cca_ord <- capscale(formula = richness ~ subset_meta$time_point* subset_meta$Population)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta$time_pointAntibiotics" = "Antibiotics",
"subset_meta$PopulationHot_dry"="Warm population",
"subset_meta$time_pointAntibiotics:subset_meta$PopulationHot_dry"="Interaction Wam population")
beta_richness_rda <-CAP_df %>%
group_by(Population, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=Population,shape = time_point)) +
scale_color_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#008080', "#d57d2c")) +
scale_fill_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#00808050', "#d57d2c50")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()
#Neutral
cca_ord <- capscale(formula = neutral ~ subset_meta$time_point* subset_meta$Population)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta$time_pointAntibiotics" = "Antibiotics",
"subset_meta$PopulationHot_dry"="Warm population",
"subset_meta$time_pointAntibiotics:subset_meta$PopulationHot_dry"="Interaction Wam population")
beta_neutral_rda <- CAP_df %>%
group_by(Population, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=Population,shape = time_point)) +
scale_color_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#008080', "#d57d2c")) +
scale_fill_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#00808050', "#d57d2c50"))+
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()
#Phylogenetic
cca_ord <- capscale(formula = phylo ~ subset_meta$time_point* subset_meta$Population)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta$time_pointAntibiotics" = "Antibiotics",
"subset_meta$PopulationHot_dry"="Warm population",
"subset_meta$time_pointAntibiotics:subset_meta$PopulationHot_dry"="Interaction Wam population")
beta_phylogenetic_rda<- CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=Population,shape = time_point)) +
scale_color_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#008080', "#d57d2c")) +
scale_fill_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#00808050', "#d57d2c50")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()ggarrange(beta_richness_rda, beta_neutral_rda, beta_phylogenetic_rda, ncol=3, nrow=1, common.legend = TRUE, legend="right")
4.3.3 Cold population
4.3.3.1 Structural zeros
acclimation_samples <- sample_metadata %>%
filter(Population == "Cold_wet" & time_point == "Acclimation") %>%
dplyr::select(sample) %>%
pull()
antibiotics_samples <- sample_metadata %>%
filter(Population == "Cold_wet" & time_point == "Antibiotics") %>%
dplyr::select(sample) %>% pull()
structural_zeros <- genome_counts_filt %>%
select(one_of(c("genome",subset_meta$sample))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
rowwise() %>% #compute for each row (genome)
mutate(all_zeros_acclimation = all(c_across(all_of(acclimation_samples)) == 0)) %>%
mutate(all_zeros_antibiotics = all(c_across(all_of(antibiotics_samples)) == 0)) %>%
mutate(average_acclimation = mean(c_across(all_of(acclimation_samples)), na.rm = TRUE)) %>%
mutate(average_antibiotics = mean(c_across(all_of(antibiotics_samples)), na.rm = TRUE)) %>%
filter(all_zeros_acclimation == TRUE || all_zeros_antibiotics==TRUE) %>%
mutate(present = case_when(
all_zeros_acclimation & !all_zeros_antibiotics ~ "antibiotics",
!all_zeros_acclimation & all_zeros_antibiotics ~ "acclimation",
!all_zeros_acclimation & !all_zeros_antibiotics ~ "None",
TRUE ~ NA_character_
)) %>%
mutate(average = ifelse(present == "acclimation", average_acclimation, average_antibiotics)) %>%
dplyr::select(genome, present, average) %>%
left_join(genome_metadata, by=join_by(genome==genome)) %>%
arrange(present,-average)Acclimation
structural_zeros %>%
filter(present=="acclimation") %>%
count(phylum, name = "sample_count") %>%
arrange(desc(sample_count)) # A tibble: 9 × 2
# Rowwise:
phylum sample_count
<chr> <int>
1 p__Bacillota_A 43
2 p__Bacteroidota 15
3 p__Bacillota 8
4 p__Pseudomonadota 7
5 p__Cyanobacteriota 3
6 p__Verrucomicrobiota 2
7 p__Bacillota_B 1
8 p__Bacillota_C 1
9 p__Fusobacteriota 1Antibiotics
# A tibble: 4 × 7
# Rowwise:
genome phylum class order family genus species
<chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 AH1_2nd_7:bin_000003 p__Pseudomonadota c__Gammaproteobacteria o__Enterobacterales f__Enterobacteriaceae g__Proteus s__Proteus cibarius
2 LI1_2nd_7:bin_000001 p__Bacillota_A c__Clostridia o__Clostridiales f__Clostridiaceae g__Sarcina s__
3 AH1_2nd_7:bin_000055 p__Bacillota c__Bacilli o__Mycoplasmatales f__Mycoplasmoidaceae g__Ureaplasma s__
4 AH1_2nd_13:bin_000025 p__Bacillota_A c__Clostridia o__Christensenellales f__UBA3700 g__ s__ Community elements differences in structural zeros
element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
element_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)# A tibble: 0 × 3
# ℹ 3 variables: trait <chr>, p_value <dbl>, p_adjust <dbl>4.3.3.2 Differential abundance analysis: composition
phylo_samples <- sample_metadata %>%
filter(Population == "Cold_wet" & time_point %in% c("Acclimation", "Antibiotics") )%>%
column_to_rownames("sample") %>%
sample_data()
phylo_genome <- genome_counts_filt %>%
filter(!genome %in% structural_zeros$genome) %>%
select(one_of(c("genome",rownames(phylo_samples)))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
column_to_rownames("genome") %>%
otu_table(., taxa_are_rows = TRUE)
phylo_taxonomy <- genome_metadata %>%
filter(genome %in% rownames(phylo_genome)) %>%
column_to_rownames("genome") %>%
dplyr::select(domain,phylum,class,order,family,genus,species) %>%
as.matrix() %>%
tax_table()
physeq_genome_filtered <- phyloseq(phylo_genome, phylo_taxonomy, phylo_samples)ancom_rand_output = ancombc2(data = physeq_genome_filtered,
assay_name = "counts",
tax_level = NULL,
fix_formula = "time_point",
rand_formula = "(1|individual)",
p_adj_method = "holm",
pseudo_sens = TRUE,
prv_cut =0,
lib_cut = 0,
s0_perc = 0.05,
group = NULL,
struc_zero = FALSE,
neg_lb = FALSE,
alpha = 0.05,
n_cl = 2,
verbose = TRUE,
global = FALSE,
pairwise = FALSE,
dunnet = FALSE,
trend = FALSE,
iter_control = list(tol = 1e-5, max_iter = 20, verbose = FALSE),
em_control = list(tol = 1e-5, max_iter = 100),
lme_control = lme4::lmerControl(),
mdfdr_control = list(fwer_ctrl_method = "holm", B = 100),
trend_control = NULL)taxonomy <- data.frame(physeq_genome_filtered@tax_table) %>%
rownames_to_column(., "taxon") %>%
mutate_at(vars(order, phylum, family, genus, species), ~ str_replace(., "[dpcofgs]__", ""))
ancombc_rand_table_mag <- ancom_rand_output$res %>%
dplyr::select(taxon, lfc_time_point2_Antibiotics, p_time_point2_Antibiotics) %>%
filter(p_time_point2_Antibiotics < 0.05) %>%
dplyr::arrange(p_time_point2_Antibiotics) %>%
merge(., taxonomy, by="taxon") %>%
mutate_at(vars(phylum, species), ~ str_replace(., "[dpcofgs]__", ""))%>%
dplyr::arrange(lfc_time_point2_Antibiotics)
colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
mutate_at(vars(phylum), ~ str_replace(., "[dpcofgs]__", "")) %>%
right_join(taxonomy, by=join_by(phylum == phylum)) %>%
dplyr::select(phylum, colors) %>%
mutate(colors = str_c(colors, "80")) %>% #add 80% alpha
unique() %>%
dplyr::arrange(phylum)
tax_table <- as.data.frame(unique(ancombc_rand_table_mag$phylum))
colnames(tax_table)[1] <- "phylum"
tax_color <- merge(tax_table, colors_alphabetic, by="phylum")%>%
dplyr::arrange(phylum) %>%
dplyr::select(colors) %>%
pull()ancombc_rand_table_mag%>%
mutate(genome=factor(taxon,levels=ancombc_rand_table_mag$taxon)) %>%
ggplot(., aes(x=lfc_time_point2_Antibiotics, y=forcats::fct_reorder(genome,lfc_time_point2_Antibiotics), fill=phylum)) + #forcats::fct_rev()
geom_col() +
scale_fill_manual(values=tax_color) +
geom_hline(yintercept=0) +
# coord_flip()+
theme(panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 14, face = "bold"),
legend.text = element_text(size = 12),
legend.title = element_text(size = 14, face = "bold"),
legend.position = "right", legend.box = "vertical")+
xlab("log2FoldChange") +
ylab("Species")+
guides(fill=guide_legend(title="Phylum"))
Phyla of the significant MAGs
ancombc_rand_table_mag%>%
filter(lfc_time_point2_Antibiotics<0) %>%
count(phylum, name = "sample_count") %>%
arrange(desc(sample_count)) phylum sample_count
1 Bacteroidota 14
2 Bacillota_A 4
3 Bacillota 2
4 Campylobacterota 1 phylum genus species
1 Campylobacterota Helicobacter_J
2 Bacillota_A
3 Bacteroidota Bacteroides
4 Bacteroidota Odoribacter
5 Bacteroidota Bacteroides
6 Bacillota Coprobacillus
7 Bacillota
8 Bacteroidota Phocaeicola
9 Bacteroidota Bacteroides
10 Bacteroidota Phocaeicola
11 Bacteroidota Odoribacter
12 Bacteroidota Bacteroides
13 Bacteroidota Bacteroides
14 Bacteroidota Parabacteroides
15 Bacteroidota
16 Bacteroidota Parabacteroides
17 Bacillota_A
18 Bacteroidota Alistipes
19 Bacillota_A
20 Bacillota_A
21 Bacteroidota Odoribacter ancombc_rand_table_mag%>%
filter(lfc_time_point2_Antibiotics<0) %>%
count(genus, name = "sample_count") %>%
arrange(desc(sample_count)) genus sample_count
1 6
2 Bacteroides 5
3 Odoribacter 3
4 Parabacteroides 2
5 Phocaeicola 2
6 Alistipes 1
7 Coprobacillus 1
8 Helicobacter_J 14.3.3.3 Differences in the functional capacity
GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts_filt$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>%
select_if(~ !is.numeric(.) || sum(.) != 0)
GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)
GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)
sample_sub <- sample_metadata %>%
filter(Population == "Cold_wet" & treatment %in% c("Acclimation", "Antibiotics"))
genome_counts_row <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>%
select(one_of(c("genome",sample_sub$sample))) %>%
column_to_rownames(., "genome")
GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_row,GIFT_db)4.3.3.3.1 Community elements differences:
MCI
GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
treatment MCI sd
<chr> <dbl> <dbl>
1 Acclimation 0.335 0.0324
2 Antibiotics 0.247 0.136 MCI <- GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.94332, p-value = 0.09304
Wilcoxon rank sum exact test
data: value by treatment
W = 190, p-value = 0.01768
alternative hypothesis: true location shift is not equal to 0element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
significant_elements <- element_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))
element_gift_t <- element_gift %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "trait")
element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "sample")%>%
left_join(., sample_metadata[c(12,13)], by = join_by(sample == sample))
difference_table <- element_gift_filt %>%
select(-sample) %>%
group_by(treatment) %>%
summarise(across(everything(), mean)) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>%
arrange(Function) %>%
mutate(Difference=Acclimation-Antibiotics)%>%
mutate(group_color = ifelse(Difference <0, "Antibiotics","Acclimation")) difference_table %>%
ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) +
geom_col() +
# geom_point(size=4) +
scale_fill_manual(values=c('#008080', '#00808050')) +
geom_hline(yintercept=0) +
coord_flip()+
theme(axis.text = element_text(size = 10),
axis.title = element_text(size = 12),
legend.position = "right",
panel.background = element_blank(),
panel.grid.major = element_line(size = 0.15, linetype = 'solid',
colour = "grey"))+
xlab("Microbial Functional Capacity") +
ylab("Mean difference")+
labs(fill="Treatment")
4.3.3.3.2 Community functions differences
MCI
GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
treatment MCI sd
<chr> <dbl> <dbl>
1 Acclimation 0.330 0.0320
2 Antibiotics 0.256 0.126 MCI <- GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.95742, p-value = 0.2331
Wilcoxon rank sum exact test
data: value by treatment
W = 188, p-value = 0.02195
alternative hypothesis: true location shift is not equal to 0function_gift <- GIFTs_functions_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
merge(., sample_metadata[c(12,13)], by="sample")unique_funct_db<- GIFT_db[c(3,4,5)] %>%
distinct(Code_function, .keep_all = TRUE)
significant_functional <- function_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
left_join(.,unique_funct_db[c(1,3)],by = join_by(trait == Code_function))| Code_function | Acclimation | Antibiotics | Function | Difference | treatment |
|---|---|---|---|---|---|
| B06 | 0.68110280 | 0.47325490 | Organic anion biosynthesis | 0.20784790 | Acclimation |
| B02 | 0.59963040 | 0.41562100 | Amino acid biosynthesis | 0.18400940 | Acclimation |
| D02 | 0.38587770 | 0.20636760 | Polysaccharide degradation | 0.17951010 | Acclimation |
| S03 | 0.27103471 | 0.09357224 | Spore | 0.17746247 | Acclimation |
| B01 | 0.84215140 | 0.69078910 | Nucleic acid biosynthesis | 0.15136230 | Acclimation |
| B07 | 0.44558700 | 0.30432910 | Vitamin biosynthesis | 0.14125790 | Acclimation |
| B08 | 0.44596150 | 0.32044710 | Aromatic compound biosynthesis | 0.12551440 | Acclimation |
| D09 | 0.25533360 | 0.13438350 | Antibiotic degradation | 0.12095010 | Acclimation |
| B04 | 0.54457570 | 0.42921780 | SCFA biosynthesis | 0.11535790 | Acclimation |
| D03 | 0.29210750 | 0.20698550 | Sugar degradation | 0.08512200 | Acclimation |
| D05 | 0.28856110 | 0.22258820 | Amino acid degradation | 0.06597290 | Acclimation |
| B10 | 0.05947806 | 0.03621687 | Antibiotic biosynthesis | 0.02326119 | Acclimation |
4.3.4 Warm population
4.3.4.1 Structural zeros
acclimation_samples <- sample_metadata %>%
filter(Population == "Hot_dry" & time_point == "Acclimation") %>%
dplyr::select(sample) %>%
pull()
antibiotics_samples <- sample_metadata %>%
filter(Population == "Hot_dry" & time_point == "Antibiotics") %>%
dplyr::select(sample) %>% pull()
structural_zeros <- genome_counts_filt %>%
select(one_of(c("genome",subset_meta$sample))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
rowwise() %>% #compute for each row (genome)
mutate(all_zeros_acclimation = all(c_across(all_of(acclimation_samples)) == 0)) %>%
mutate(all_zeros_antibiotics = all(c_across(all_of(antibiotics_samples)) == 0)) %>%
mutate(average_acclimation = mean(c_across(all_of(acclimation_samples)), na.rm = TRUE)) %>%
mutate(average_antibiotics = mean(c_across(all_of(antibiotics_samples)), na.rm = TRUE)) %>%
filter(all_zeros_acclimation == TRUE || all_zeros_antibiotics==TRUE) %>%
mutate(present = case_when(
all_zeros_acclimation & !all_zeros_antibiotics ~ "antibiotics",
!all_zeros_acclimation & all_zeros_antibiotics ~ "acclimation",
!all_zeros_acclimation & !all_zeros_antibiotics ~ "None",
TRUE ~ NA_character_
)) %>%
mutate(average = ifelse(present == "acclimation", average_acclimation, average_antibiotics)) %>%
dplyr::select(genome, present, average) %>%
left_join(genome_metadata, by=join_by(genome==genome)) %>%
arrange(present,-average)Acclimation
structural_zeros %>%
filter(present=="acclimation") %>%
count(phylum, name = "sample_count") %>%
arrange(desc(sample_count)) # A tibble: 12 × 2
# Rowwise:
phylum sample_count
<chr> <int>
1 p__Bacillota_A 63
2 p__Bacteroidota 32
3 p__Bacillota 16
4 p__Cyanobacteriota 9
5 p__Pseudomonadota 7
6 p__Bacillota_B 3
7 p__Desulfobacterota 3
8 p__Bacillota_C 2
9 p__Verrucomicrobiota 2
10 p__Actinomycetota 1
11 p__Campylobacterota 1
12 p__Elusimicrobiota 1Antibiotics
# A tibble: 7 × 7
# Rowwise:
genome phylum class order family genus species
<chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 AH1_2nd_7:bin_000003 p__Pseudomonadota c__Gammaproteobacteria o__Enterobacterales f__Enterobacteriaceae g__Proteus s__Proteus cibarius
2 LI1_2nd_8:bin_000030 p__Actinomycetota c__Actinomycetia o__Mycobacteriales f__Mycobacteriaceae g__Corynebacterium s__
3 LI1_2nd_8:bin_000077 p__Bacteroidota c__Bacteroidia o__Bacteroidales f__Tannerellaceae g__Parabacteroides s__Parabacteroides sp003480915
4 LI1_2nd_5:bin_000032 p__Bacillota_A c__Clostridia o__Christensenellales f__UBA1242 g__Caccovivens s__
5 AH1_2nd_18:bin_000013 p__Bacteroidota c__Bacteroidia o__Bacteroidales f__Tannerellaceae g__Parabacteroides s__
6 LI1_2nd_5:bin_000069 p__Bacillota c__Bacilli o__RF39 f__UBA660 g__ s__
7 AH1_2nd_20:bin_000061 p__Bacillota c__Bacilli o__Lactobacillales f__Enterococcaceae g__Enterococcus s__ Community elements differences in structural zeros
element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
element_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)# A tibble: 65 × 3
trait p_value p_adjust
<chr> <dbl> <dbl>
1 B0101 0.00280 0.0184
2 B0103 0.0196 0.0437
3 B0105 0.00559 0.0184
4 B0106 0.00559 0.0184
5 B0205 0.00280 0.0184
6 B0207 0.0112 0.0280
7 B0212 0.00559 0.0184
8 B0214 0.00662 0.0184
9 B0216 0.00559 0.0184
10 B0307 0.0196 0.0437
# ℹ 55 more rows4.3.4.2 Differential abundance analysis: composition
phylo_samples <- sample_metadata %>%
filter(Population == "Hot_dry" & time_point %in% c("Acclimation", "Antibiotics") )%>%
column_to_rownames("sample") %>%
sample_data()
phylo_genome <- genome_counts_filt %>%
filter(!genome %in% structural_zeros$genome) %>%
select(one_of(c("genome",rownames(phylo_samples)))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
column_to_rownames("genome") %>%
otu_table(., taxa_are_rows = TRUE)
phylo_taxonomy <- genome_metadata %>%
filter(genome %in% rownames(phylo_genome)) %>%
column_to_rownames("genome") %>%
dplyr::select(domain,phylum,class,order,family,genus,species) %>%
as.matrix() %>%
tax_table()
physeq_genome_filtered <- phyloseq(phylo_genome, phylo_taxonomy, phylo_samples)ancom_rand_output = ancombc2(data = physeq_genome_filtered,
assay_name = "counts",
tax_level = NULL,
fix_formula = "time_point",
rand_formula = "(1|individual)",
p_adj_method = "holm",
pseudo_sens = TRUE,
prv_cut =0,
lib_cut = 0,
s0_perc = 0.05,
group = NULL,
struc_zero = FALSE,
neg_lb = FALSE,
alpha = 0.05,
n_cl = 2,
verbose = TRUE,
global = FALSE,
pairwise = FALSE,
dunnet = FALSE,
trend = FALSE,
iter_control = list(tol = 1e-5, max_iter = 20, verbose = FALSE),
em_control = list(tol = 1e-5, max_iter = 100),
lme_control = lme4::lmerControl(),
mdfdr_control = list(fwer_ctrl_method = "holm", B = 100),
trend_control = NULL)taxonomy <- data.frame(physeq_genome_filtered@tax_table) %>%
rownames_to_column(., "taxon") %>%
mutate_at(vars(order, phylum, family, genus, species), ~ str_replace(., "[dpcofgs]__", ""))
ancombc_rand_table_mag <- ancom_rand_output$res %>%
dplyr::select(taxon, lfc_time_point2_Antibiotics, p_time_point2_Antibiotics) %>%
filter(p_time_point2_Antibiotics < 0.05) %>%
dplyr::arrange(p_time_point2_Antibiotics) %>%
merge(., taxonomy, by="taxon") %>%
mutate_at(vars(phylum, species), ~ str_replace(., "[dpcofgs]__", ""))%>%
dplyr::arrange(lfc_time_point2_Antibiotics)
colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
mutate_at(vars(phylum), ~ str_replace(., "[dpcofgs]__", "")) %>%
right_join(taxonomy, by=join_by(phylum == phylum)) %>%
dplyr::select(phylum, colors) %>%
mutate(colors = str_c(colors, "80")) %>% #add 80% alpha
unique() %>%
dplyr::arrange(phylum)
tax_table <- as.data.frame(unique(ancombc_rand_table_mag$phylum))
colnames(tax_table)[1] <- "phylum"
tax_color <- merge(tax_table, colors_alphabetic, by="phylum")%>%
dplyr::arrange(phylum) %>%
dplyr::select(colors) %>%
pull()ancombc_rand_table_mag%>%
mutate(genome=factor(taxon,levels=ancombc_rand_table_mag$taxon)) %>%
ggplot(., aes(x=lfc_time_point2_Antibiotics, y=forcats::fct_reorder(genome,lfc_time_point2_Antibiotics), fill=phylum)) + #forcats::fct_rev()
geom_col() +
scale_fill_manual(values=tax_color) +
geom_hline(yintercept=0) +
# coord_flip()+
theme(panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 14, face = "bold"),
legend.text = element_text(size = 12),
legend.title = element_text(size = 14, face = "bold"),
legend.position = "right", legend.box = "vertical")+
xlab("log2FoldChange") +
ylab("Species")+
guides(fill=guide_legend(title="Phylum"))Phyla of the significant MAGs
ancombc_rand_table_mag%>%
filter(lfc_time_point2_Antibiotics<0) %>%
count(phylum, name = "sample_count") %>%
arrange(desc(sample_count))
ancombc_rand_table_mag%>%
filter(lfc_time_point2_Antibiotics<0) %>%
select(phylum, genus, species)
ancombc_rand_table_mag%>%
filter(lfc_time_point2_Antibiotics<0) %>%
count(genus, name = "sample_count") %>%
arrange(desc(sample_count))4.3.4.3 Differences in the functional capacity
GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts_filt$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>%
select_if(~ !is.numeric(.) || sum(.) != 0)
GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)
GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)
sample_sub <- sample_metadata %>%
filter(Population == "Cold_wet" & treatment %in% c("Acclimation", "Antibiotics"))
genome_counts_row <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>%
select(one_of(c("genome",sample_sub$sample))) %>%
column_to_rownames(., "genome")
GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_row,GIFT_db)4.3.4.3.1 Community elements differences:
MCI
GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
treatment MCI sd
<chr> <dbl> <dbl>
1 Acclimation 0.335 0.0324
2 Antibiotics 0.247 0.136 MCI <- GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.94332, p-value = 0.09304
Wilcoxon rank sum exact test
data: value by treatment
W = 190, p-value = 0.01768
alternative hypothesis: true location shift is not equal to 0element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
significant_elements <- element_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))
element_gift_t <- element_gift %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "trait")
element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "sample")%>%
left_join(., sample_metadata[c(12,13)], by = join_by(sample == sample))
difference_table <- element_gift_filt %>%
select(-sample) %>%
group_by(treatment) %>%
summarise(across(everything(), mean)) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>%
arrange(Function) %>%
mutate(Difference=Acclimation-Antibiotics)%>%
mutate(group_color = ifelse(Difference <0, "Antibiotics","Acclimation")) difference_table %>%
ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) +
geom_col() +
# geom_point(size=4) +
scale_fill_manual(values=c("#d57d2c50","#d57d2c")) +
geom_hline(yintercept=0) +
coord_flip()+
theme(axis.text = element_text(size = 10),
axis.title = element_text(size = 12),
legend.position = "right",
panel.background = element_blank(),
panel.grid.major = element_line(size = 0.15, linetype = 'solid',
colour = "grey"))+
xlab("Microbial Functional Capacity") +
ylab("Mean difference")+
labs(fill="Treatment")
4.3.4.3.2 Community functions differences
MCI
GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
treatment MCI sd
<chr> <dbl> <dbl>
1 Acclimation 0.330 0.0320
2 Antibiotics 0.256 0.126 MCI <- GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.95742, p-value = 0.2331
Wilcoxon rank sum exact test
data: value by treatment
W = 188, p-value = 0.02195
alternative hypothesis: true location shift is not equal to 0function_gift <- GIFTs_functions_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
merge(., sample_metadata[c(12,13)], by="sample")unique_funct_db<- GIFT_db[c(3,4,5)] %>%
distinct(Code_function, .keep_all = TRUE)
significant_functional <- function_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
left_join(.,unique_funct_db[c(1,3)],by = join_by(trait == Code_function))| Code_function | Acclimation | Antibiotics | Function | Difference | treatment |
|---|---|---|---|---|---|
| B06 | 0.68110280 | 0.47325490 | Organic anion biosynthesis | 0.20784790 | Acclimation |
| B02 | 0.59963040 | 0.41562100 | Amino acid biosynthesis | 0.18400940 | Acclimation |
| D02 | 0.38587770 | 0.20636760 | Polysaccharide degradation | 0.17951010 | Acclimation |
| S03 | 0.27103471 | 0.09357224 | Spore | 0.17746247 | Acclimation |
| B01 | 0.84215140 | 0.69078910 | Nucleic acid biosynthesis | 0.15136230 | Acclimation |
| B07 | 0.44558700 | 0.30432910 | Vitamin biosynthesis | 0.14125790 | Acclimation |
| B08 | 0.44596150 | 0.32044710 | Aromatic compound biosynthesis | 0.12551440 | Acclimation |
| D09 | 0.25533360 | 0.13438350 | Antibiotic degradation | 0.12095010 | Acclimation |
| B04 | 0.54457570 | 0.42921780 | SCFA biosynthesis | 0.11535790 | Acclimation |
| D03 | 0.29210750 | 0.20698550 | Sugar degradation | 0.08512200 | Acclimation |
| D05 | 0.28856110 | 0.22258820 | Amino acid degradation | 0.06597290 | Acclimation |
| B10 | 0.05947806 | 0.03621687 | Antibiotic biosynthesis | 0.02326119 | Acclimation |
4.4 Does antibiotic administration remove the differences found in the two populations?
4.4.1 Shapiro test
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point == "Antibiotics" ) %>%
filter(metric=="richness") %>%
summarize(shapiro_p_value = shapiro.test(value)$p.value) %>%
pull(shapiro_p_value)[1] 0.001165318alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point == "Antibiotics" ) %>%
filter(metric=="neutral") %>%
summarize(shapiro_p_value = shapiro.test(value)$p.value) %>%
pull(shapiro_p_value)[1] 0.0003462674alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point == "Antibiotics" ) %>%
filter(metric=="phylogenetic") %>%
summarize(shapiro_p_value = shapiro.test(value)$p.value) %>%
pull(shapiro_p_value)[1] 0.066284424.4.2 Alpha diversity
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == sample))%>%
filter(time_point == "Antibiotics" )alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
mutate(metric=factor(metric,levels=c("richness", "neutral"))) %>%
filter(metric!="phylogenetic") %>%
filter(time_point == "Antibiotics" ) %>%
ggplot(aes(y = value, x = Population, color=Population, fill=Population)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#008080', "#d57d2c")) +
scale_fill_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#00808050', "#d57d2c50")) +
stat_compare_means(method = "t.test", show.legend = F, size = 3, label.x = c(1.5))+
facet_grid( ~ metric, scales = "free_y")+
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
mutate(metric=factor(metric,levels=c("richness","neutral"))) %>%
filter(metric!="phylogenetic") %>%
filter(time_point == "Antibiotics" ) %>%
ggplot(aes(y = value, x = Population, color=Population, fill=Population)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#008080', "#d57d2c")) +
scale_fill_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#00808050', "#d57d2c50")) +
stat_compare_means(method = "wilcox.test", show.legend = F, size = 3, label.x = c(1.5))+
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")
4.4.3 Beta diversity
samples_to_keep <- sample_metadata %>%
filter(time_point == "Antibiotics") %>%
select(Tube_code) %>%
pull()
subset_meta <- sample_metadata %>%
filter(time_point == "Antibiotics")Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$Population) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.015377 0.0153774 6.8934 999 0.021 *
Residuals 21 0.046845 0.0022307
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(richness ~ Population,
data = subset_meta %>% arrange(match(Tube_code,labels(richness))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 1.361743 | 0.1533845 | 3.80465 | 0.001 |
| Residual | 21 | 7.516224 | 0.8466155 | NA | NA |
| Total | 22 | 8.877966 | 1.0000000 | NA | NA |
Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$Population) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.030649 0.0306488 3.8694 999 0.06 .
Residuals 21 0.166339 0.0079209
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(neutral ~ Population,
data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 1.787698 | 0.208772 | 5.541022 | 0.001 |
| Residual | 21 | 6.775221 | 0.791228 | NA | NA |
| Total | 22 | 8.562919 | 1.000000 | NA | NA |
Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$Population) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.012165 0.012165 0.9979 999 0.33
Residuals 21 0.256012 0.012191 adonis2(phylo ~ Population,
data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.8970751 | 0.1890846 | 4.896659 | 0.001 |
| Residual | 21 | 3.8472307 | 0.8109154 | NA | NA |
| Total | 22 | 4.7443057 | 1.0000000 | NA | NA |
4.5 Are the microbial communities similar in both donor samples?
4.5.1 Alpha diversity
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == sample))%>%
filter(type =="Hot_control" & time_point %in% c( "Transplant1", "Transplant2"))samples_to_keep <- sample_metadata %>%
filter(type =="Hot_control" & time_point %in% c( "Transplant1", "Transplant2"))%>%
select(sample) %>%
pull()
subset_meta <- sample_metadata %>%
filter(type =="Hot_control" & time_point %in% c( "Transplant1", "Transplant2"))alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic"))) %>%
filter(type =="Hot_control" & time_point %in% c( "Transplant1", "Transplant2")) %>%
ggplot(aes(y = value, x = time_point, color=time_point, fill=time_point)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="time_point",
breaks=c("Transplant1","Transplant2"),
values=c('#d5992c', "#d5b52c")) +
scale_fill_manual(name="time_point",
breaks=c("Transplant1","Transplant2"),
values=c('#d5992c50', "#d5b52c50")) +
facet_wrap(. ~ metric, scales = "free") +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")
Richness
Modelq0GLMMNB <- glmer.nb(richness ~ time_point+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Negative Binomial(17.9668) ( log )
Formula: richness ~ time_point + (1 | individual)
Data: alpha_div_meta
AIC BIC logLik deviance df.resid
150.4 153.3 -71.2 142.4 11
Scaled residuals:
Min 1Q Median 3Q Max
-2.02956 -0.52933 0.03513 0.56378 1.56130
Random effects:
Groups Name Variance Std.Dev.
individual (Intercept) 0.009989 0.09995
Number of obs: 15, groups: individual, 8
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.60703 0.09905 46.514 <2e-16 ***
time_pointTransplant2 0.05482 0.13299 0.412 0.68
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
tm_pntTrns2 -0.624Neutral
Model_neutral <- nlme::lme(fixed = neutral ~ time_point, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_neutral)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
115.9183 118.1781 -53.95914
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 9.533528 10.15744
Fixed effects: neutral ~ time_point
Value Std.Error DF t-value p-value
(Intercept) 43.94053 4.925212 7 8.921551 0.0000
time_pointTransplant2 4.31623 5.338413 6 0.808523 0.4497
Correlation:
(Intr)
time_pointTransplant2 -0.491
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.57714267 -0.56764018 0.05406399 0.55332481 1.12750552
Number of Observations: 15
Number of Groups: 8 Phylogenetic
Model_phylo <- nlme::lme(fixed = phylogenetic ~ time_point, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_phylo)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
41.16537 43.42517 -16.58269
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 1.16159 0.3042944
Fixed effects: phylogenetic ~ time_point
Value Std.Error DF t-value p-value
(Intercept) 5.785003 0.4245418 7 13.626462 0.0000
time_pointTransplant2 -0.018098 0.1623255 6 -0.111494 0.9149
Correlation:
(Intr)
time_pointTransplant2 -0.168
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.2428123 -0.3756145 -0.0502414 0.4569146 1.3175579
Number of Observations: 15
Number of Groups: 8 4.5.2 Beta diversity
Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.00014 0.0001396 0.0173 999 0.907
Residuals 13 0.10481 0.0080621 adonis2(richness ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(richness))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(richness))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.06889268 | 0.03010712 | 0.4035421 | 0.9375 |
| Residual | 13 | 2.21935895 | 0.96989288 | NA | NA |
| Total | 14 | 2.28825162 | 1.00000000 | NA | NA |
Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.000655 0.0006546 0.0514 999 0.835
Residuals 13 0.165409 0.0127237 adonis2(neutral ~ treatment,
data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(neutral))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.09294671 | 0.03882177 | 0.5250671 | 0.7265625 |
| Residual | 13 | 2.30124351 | 0.96117823 | NA | NA |
| Total | 14 | 2.39419022 | 1.00000000 | NA | NA |
Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.003691 0.0036912 0.3071 999 0.674
Residuals 13 0.156266 0.0120205 adonis2(phylo ~ treatment,
data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(phylo))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.009773632 | 0.01921019 | 0.2546239 | 0.7734375 |
| Residual | 13 | 0.498999565 | 0.98078981 | NA | NA |
| Total | 14 | 0.508773196 | 1.00000000 | NA | NA |
4.6 Does the donor sample maintain the microbial community found in acclimation?
sample_metadata <- sample_metadata %>%
mutate(treatment=case_when(
treatment %in% c("Transplant1", "Transplant2") ~ "Transplant",
TRUE ~ treatment
))
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == sample))%>%
filter(type == "Hot_control" & treatment %in% c("Acclimation","Transplant"))sample_metadata <- sample_metadata %>%
mutate(treatment=case_when(
treatment %in% c("Transplant1", "Transplant2") ~ "Transplant",
TRUE ~ treatment
))
samples_to_keep <- sample_metadata %>%
filter(type == "Hot_control" & treatment %in% c("Acclimation","Transplant")) %>%
select(sample) %>%
pull()
subset_meta <- sample_metadata %>%
filter(type == "Hot_control" & treatment %in% c("Acclimation","Transplant"))4.6.1 Alpha diversity
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic"))) %>%
filter(type == "Hot_control" & treatment %in% c("Acclimation","Transplant")) %>%
ggplot(aes(y = value, x = treatment, color=treatment, fill=treatment)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="treatment",
breaks=c("Acclimation","Transplant"),
values=c("#d57d2c", "#d5b52c")) +
scale_fill_manual(name="treatment",
breaks=c("Acclimation","Transplant"),
values=c("#d57d2c50", "#d5b52c50")) +
facet_wrap(. ~ metric, scales = "free") +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")
Richness
Modelq0GLMMNB <- glmer.nb(richness ~ treatment+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Negative Binomial(19.6072) ( log )
Formula: richness ~ treatment + (1 | individual)
Data: alpha_div_meta
AIC BIC logLik deviance df.resid
228.0 232.7 -110.0 220.0 20
Scaled residuals:
Min 1Q Median 3Q Max
-2.37872 -0.53180 0.06467 0.46904 1.76837
Random effects:
Groups Name Variance Std.Dev.
individual (Intercept) 0 0
Number of obs: 24, groups: individual, 9
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.4569 0.0834 53.441 <2e-16 ***
treatmentTransplant 0.1855 0.1049 1.769 0.0769 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
trtmntTrnsp -0.795
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')Neutral
Model_neutral <- nlme::lme(fixed = neutral ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_neutral)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
185.3367 189.7009 -88.66837
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 0.001909967 12.18199
Fixed effects: neutral ~ treatment
Value Std.Error DF t-value p-value
(Intercept) 44.09058 4.060663 14 10.857976 0.0000
treatmentTransplant 1.80350 5.136377 14 0.351122 0.7307
Correlation:
(Intr)
treatmentTransplant -0.791
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.66610957 -0.48314823 -0.03902931 0.62479437 2.02597977
Number of Observations: 24
Number of Groups: 9 Phylogenetic
Model_phylo <- nlme::lme(fixed = phylogenetic ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_phylo)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
80.6646 85.02877 -36.3323
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 0.7197786 0.9593315
Fixed effects: phylogenetic ~ treatment
Value Std.Error DF t-value p-value
(Intercept) 6.463020 0.3997775 14 16.166544 0.0000
treatmentTransplant -0.577492 0.4112749 14 -1.404151 0.1821
Correlation:
(Intr)
treatmentTransplant -0.622
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.329685717 -0.568101055 -0.001989187 0.552370170 1.529548991
Number of Observations: 24
Number of Groups: 9 4.6.2 Beta diversity
Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.00836 0.0083600 1.4409 999 0.273
Residuals 22 0.12764 0.0058019 adonis2(richness ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(richness))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(richness))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.2657351 | 0.06396199 | 1.503319 | 0.111 |
| Residual | 22 | 3.8888430 | 0.93603801 | NA | NA |
| Total | 23 | 4.1545781 | 1.00000000 | NA | NA |
Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.000661 0.0006614 0.0736 999 0.798
Residuals 22 0.197804 0.0089911 adonis2(neutral ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(neutral))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(neutral))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.3164816 | 0.07596593 | 1.808646 | 0.078 |
| Residual | 22 | 3.8496178 | 0.92403407 | NA | NA |
| Total | 23 | 4.1660995 | 1.00000000 | NA | NA |
Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.00204 0.0020405 0.2622 999 0.683
Residuals 22 0.17118 0.0077807 adonis2(phylo ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(phylo))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(phylo))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.0412056 | 0.056244 | 1.31111 | 0.235 |
| Residual | 22 | 0.6914166 | 0.943756 | NA | NA |
| Total | 23 | 0.7326222 | 1.000000 | NA | NA |
p0<-richness %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
scale_color_manual(name="treatment",
breaks=c("Acclimation","Transplant"),
values=c("#d57d2c", "#d5b52c")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)
p1<-neutral %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
scale_color_manual(name="treatment",
breaks=c("Acclimation","Transplant"),
values=c("#d57d2c", "#d5b52c")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)
p2<-phylo %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
scale_color_manual(name="treatment",
breaks=c("Acclimation","Transplant"),
values=c("#d57d2c", "#d5b52c")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)
4.7 Is the donor sample microbiota different to recipient microbiota?
samples_to_keep <- sample_metadata %>%
filter(type == "Treatment" & treatment %in% c("Acclimation", "Transplant")) %>%
select(sample) %>%
pull()
subset_meta <- sample_metadata %>%
filter(type == "Treatment" & treatment %in% c("Acclimation", "Transplant"))4.7.1 Alpha diversity
sample_metadata <- sample_metadata %>%
mutate(treatment=case_when(
treatment %in% c("Transplant1", "Transplant2") ~ "Transplant",
TRUE ~ treatment
))
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == sample))%>%
filter(type == "Treatment" & treatment %in% c("Acclimation","Transplant"))alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic"))) %>%
filter(type == "Treatment" & treatment %in% c("Acclimation", "Transplant")) %>%
ggplot(aes(y = value, x = treatment, color=treatment, fill=treatment)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="treatment",
breaks=c("Acclimation","Transplant"),
values=c('#008080', "#d5b52c")) +
scale_fill_manual(name="treatment",
breaks=c("Acclimation","Transplant"),
values=c('#00808050', "#d5b52c50")) +
facet_wrap(. ~ metric, scales = "free") +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")
Richness
Modelq0GLMMNB <- glmer.nb(richness ~ treatment+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Negative Binomial(6.4835) ( log )
Formula: richness ~ treatment + (1 | individual)
Data: alpha_div_meta
AIC BIC logLik deviance df.resid
219.3 223.6 -105.6 211.3 18
Scaled residuals:
Min 1Q Median 3Q Max
-1.9573 -0.5447 0.1304 0.5773 1.2167
Random effects:
Groups Name Variance Std.Dev.
individual (Intercept) 5.302e-13 7.281e-07
Number of obs: 22, groups: individual, 9
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.8067 0.1400 27.19 < 2e-16 ***
treatmentTransplant 0.8399 0.1795 4.68 2.87e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
trtmntTrnsp -0.780
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')Neutral
Model_neutral <-nlme:: lme(fixed = neutral ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_neutral)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
169.583 173.5659 -80.79148
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 4.713488 11.40939
Fixed effects: neutral ~ treatment
Value Std.Error DF t-value p-value
(Intercept) 17.31015 4.114893 12 4.206708 0.0012
treatmentTransplant 28.70744 5.016259 12 5.722879 0.0001
Correlation:
(Intr)
treatmentTransplant -0.701
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.4861810 -0.5958348 0.1169620 0.5490436 1.7889456
Number of Observations: 22
Number of Groups: 9 Phylogenetic
Model_phylo <- nlme::lme(fixed = phylogenetic ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_phylo)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
80.93175 84.91468 -36.46587
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 0.8689563 1.118803
Fixed effects: phylogenetic ~ treatment
Value Std.Error DF t-value p-value
(Intercept) 5.488256 0.4722057 12 11.622595 0.00
treatmentTransplant 0.226402 0.5020134 12 0.450988 0.66
Correlation:
(Intr)
treatmentTransplant -0.587
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.42844864 -0.44112019 -0.05260838 0.51238319 1.55109109
Number of Observations: 22
Number of Groups: 9 4.7.2 Beta diversity
Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.12795 0.127946 13.179 999 0.001 ***
Residuals 20 0.19416 0.009708
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(richness ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(richness))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(richness))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 1.777815 | 0.2760196 | 7.625059 | 0.001 |
| Residual | 20 | 4.663086 | 0.7239804 | NA | NA |
| Total | 21 | 6.440902 | 1.0000000 | NA | NA |
Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.071502 0.071502 5.4967 999 0.029 *
Residuals 20 0.260166 0.013008
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(neutral ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(neutral))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(neutral))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 2.112572 | 0.3210813 | 9.458609 | 0.001 |
| Residual | 20 | 4.466983 | 0.6789187 | NA | NA |
| Total | 21 | 6.579556 | 1.0000000 | NA | NA |
Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.04731 0.047307 2.6157 999 0.117
Residuals 20 0.36172 0.018086 adonis2(phylo ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(phylo))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(phylo))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.3778248 | 0.239656 | 6.303884 | 0.001 |
| Residual | 20 | 1.1987049 | 0.760344 | NA | NA |
| Total | 21 | 1.5765298 | 1.000000 | NA | NA |
p0<-richness %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
scale_color_manual(name="treatment",
breaks=c("Acclimation","Transplant"),
values=c('#008080', "#d5b52c")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)
p1<-neutral %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
scale_color_manual(name="treatment",
breaks=c("Acclimation","Transplant"),
values=c('#008080', "#d5b52c")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)
p2<-phylo %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
scale_color_manual(name="treatment",
breaks=c("Acclimation","Transplant"),
values=c('#008080', "#d5b52c")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)
4.8 Does FMT change the microbial community over time?
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == sample))%>%
filter(type == "Treatment" & treatment %in% c("FMT1","FMT2") ) 4.8.1 Alpha diversity
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
mutate(metric=factor(metric,levels=c("richness","neutral", "phylogenetic"))) %>%
filter(type=="Treatment" & treatment %in% c("FMT1", "FMT2")) %>%
ggplot(aes(y = value, x = treatment, color=treatment, fill=treatment)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="treatment",
values=c("#76b183", "#2b4b32")) +
scale_fill_manual(name="treatment",
values=c("#76b18350", "#2b4b3250")) +
facet_wrap(. ~ metric, scales = "free") +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")
Richness
Modelq0GLMMNB <- glmer.nb(richness ~ treatment+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Negative Binomial(4.3876) ( log )
Formula: richness ~ treatment + (1 | individual)
Data: alpha_div_meta
AIC BIC logLik deviance df.resid
171.0 174.3 -81.5 163.0 13
Scaled residuals:
Min 1Q Median 3Q Max
-1.73495 -0.71265 -0.07086 0.40744 1.88240
Random effects:
Groups Name Variance Std.Dev.
individual (Intercept) 1.073e-11 3.275e-06
Number of obs: 17, groups: individual, 9
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.9343 0.1759 22.369 <2e-16 ***
treatmentFMT2 0.4052 0.2402 1.687 0.0917 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
tretmntFMT2 -0.732
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')Neutral
Model_neutral <- nlme::lme(fixed = neutral ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_neutral)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
128.7946 131.6268 -60.3973
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 3.497472 11.25384
Fixed effects: neutral ~ treatment
Value Std.Error DF t-value p-value
(Intercept) 24.65947 4.164756 8 5.920988 0.0004
treatmentFMT2 14.87494 5.482531 7 2.713151 0.0301
Correlation:
(Intr)
treatmentFMT2 -0.7
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.76462324 -0.55701822 -0.04763166 0.55267588 1.30333436
Number of Observations: 17
Number of Groups: 9 Phylogenetic
Model_phylo <- nlme::lme(fixed = phylogenetic ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_phylo)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
56.05281 58.88501 -24.0264
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 0.6954001 0.8350046
Fixed effects: phylogenetic ~ treatment
Value Std.Error DF t-value p-value
(Intercept) 4.154090 0.3805932 8 10.914778 0.0000
treatmentFMT2 0.912624 0.4105974 7 2.222673 0.0616
Correlation:
(Intr)
treatmentFMT2 -0.583
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.1361565 -0.5779234 -0.1465750 0.4815529 2.3352020
Number of Observations: 17
Number of Groups: 9 4.8.2 Beta diversity
samples_to_keep <- sample_metadata %>%
filter(type == "Treatment" & treatment %in% c("FMT1", "FMT2")) %>%
select(sample) %>%
pull()
subset_meta <- sample_metadata %>%
filter(type == "Treatment" & treatment %in% c("FMT1", "FMT2"))Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.017610 0.017610 2.8225 999 0.115
Residuals 15 0.093585 0.006239 adonis2(richness ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(richness))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(richness))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.3560084 | 0.07968734 | 1.298809 | 0.00390625 |
| Residual | 15 | 4.1115570 | 0.92031266 | NA | NA |
| Total | 16 | 4.4675655 | 1.00000000 | NA | NA |
Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.009828 0.0098278 0.7113 999 0.442
Residuals 15 0.207263 0.0138175 adonis2(neutral ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(neutral))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(neutral))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.3149954 | 0.08110541 | 1.323962 | 0.00390625 |
| Residual | 15 | 3.5687823 | 0.91889459 | NA | NA |
| Total | 16 | 3.8837776 | 1.00000000 | NA | NA |
Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.010955 0.010955 0.6602 999 0.52
Residuals 15 0.248892 0.016593 adonis2(phylo ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(phylo))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(phylo))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.06391536 | 0.09449338 | 1.565312 | 0.0703125 |
| Residual | 15 | 0.61248501 | 0.90550662 | NA | NA |
| Total | 16 | 0.67640037 | 1.00000000 | NA | NA |
p0<-richness %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
scale_color_manual(name="treatment",
breaks=c("FMT1", "FMT2"),
values=c("#76b183", "#2b4b32")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)
p1<-neutral %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
scale_color_manual(name="treatment",
breaks=c("FMT1", "FMT2"),
values=c("#76b183", "#2b4b32")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)
p2<-phylo %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
scale_color_manual(name="treatment",
breaks=c("FMT1", "FMT2"),
values=c("#76b183", "#2b4b32")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)
4.9 Do FMT change the microbial community compared to antibiotics and acclimation?
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == sample))%>%
filter(type == "Treatment" & treatment %in% c("Antibiotics", "Acclimation","FMT1") ) 4.9.1 Alpha diversity
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
mutate(metric=factor(metric,levels=c("richness","neutral", "phylogenetic"))) %>%
filter(type=="Treatment" & treatment %in% c( "Antibiotics","Acclimation", "FMT1")) %>%
ggplot(aes(y = value, x = treatment, color=treatment, fill=treatment)) +
geom_jitter(width = 0.2, show.legend = FALSE) +
geom_boxplot(width = 0.5, alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
scale_color_manual(name="Treatment",
breaks=c("Antibiotics","Acclimation", "FMT1"),
values=c("#003a45",'#008080', "#76b183")) +
scale_fill_manual(name="Treatment",
breaks=c("Antibiotics","Acclimation", "FMT1"),
values=c("#003a4550",'#00808050',"#76b18350")) +
facet_wrap(. ~ metric, scales = "free") +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")
Richness: Problems to calculate
Modelq0GLMMNB <- glmer.nb(richness ~ treatment+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)Neutral
Model_neutral <- lme(fixed = neutral ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_neutral)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
173.0404 178.263 -81.5202
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 4.263247 9.322271
Fixed effects: neutral ~ treatment
Value Std.Error DF t-value p-value
(Intercept) 17.310151 3.416951 13 5.065963 0.0002
treatmentAntibiotics -8.868175 4.744329 13 -1.869216 0.0843
treatmentFMT1 7.289358 4.552796 13 1.601073 0.1334
Correlation:
(Intr) trtmnA
treatmentAntibiotics -0.596
treatmentFMT1 -0.621 0.457
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.93185401 -0.56384573 -0.04015247 0.47657219 1.55593253
Number of Observations: 24
Number of Groups: 9 Phylogenetic
Model_phylo <- lme(fixed = phylogenetic ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_phylo)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
92.7946 98.01721 -41.3973
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 0.5507824 1.405346
Fixed effects: phylogenetic ~ treatment
Value Std.Error DF t-value p-value
(Intercept) 5.488256 0.5031412 13 10.907984 0.0000
treatmentAntibiotics -1.198294 0.7137109 13 -1.678962 0.1170
treatmentFMT1 -1.351708 0.6855445 13 -1.971729 0.0703
Correlation:
(Intr) trtmnA
treatmentAntibiotics -0.611
treatmentFMT1 -0.636 0.456
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.00120146 -0.44426020 0.00785757 0.38719755 1.97397308
Number of Observations: 24
Number of Groups: 9 4.9.2 Beta diversity
samples_to_keep <- sample_metadata %>%
filter(type == "Treatment" & treatment %in% c("Acclimation", "Antibiotics", "FMT1")) %>%
select(sample) %>%
pull()
subset_meta <- sample_metadata %>%
filter(type == "Treatment" & treatment %in% c("Acclimation", "Antibiotics", "FMT1"))Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.01582 0.0079101 1.05 999 0.397
Residuals 21 0.15821 0.0075336 adonis2(richness ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(richness))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(richness))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 2 | 1.840841 | 0.2018508 | 2.655435 | 0.001 |
| Residual | 21 | 7.278967 | 0.7981492 | NA | NA |
| Total | 23 | 9.119808 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Acclimation vs Antibiotics 1 0.8287871 2.293722 0.1407734 0.003 0.009 *
2 Acclimation vs FMT1 1 0.9572484 2.939042 0.1638350 0.001 0.003 *
3 Antibiotics vs FMT1 1 0.9764237 2.751191 0.1746656 0.006 0.018 .Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.01796 0.0089778 0.5959 999 0.537
Residuals 21 0.31640 0.0150667 adonis2(neutral ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(neutral))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(neutral))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 2 | 2.357142 | 0.2699241 | 3.882065 | 0.001 |
| Residual | 21 | 6.375470 | 0.7300759 | NA | NA |
| Total | 23 | 8.732613 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Acclimation vs Antibiotics 1 0.8814932 2.747044 0.1640316 0.001 0.003 *
2 Acclimation vs FMT1 1 1.3133104 4.634354 0.2360329 0.001 0.003 *
3 Antibiotics vs FMT1 1 1.3427497 4.355528 0.2509591 0.001 0.003 *Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.15410 0.077052 2.6434 999 0.102
Residuals 21 0.61214 0.029149 adonis2(phylo ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(phylo))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(phylo))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 2 | 1.327076 | 0.3647564 | 6.029092 | 0.001 |
| Residual | 21 | 2.311178 | 0.6352436 | NA | NA |
| Total | 23 | 3.638254 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Acclimation vs Antibiotics 1 0.6885926 5.124832 0.2679674 0.004 0.012 .
2 Acclimation vs FMT1 1 0.2548027 3.292748 0.1800029 0.040 0.120
3 Antibiotics vs FMT1 1 1.1000473 9.048069 0.4103792 0.004 0.012 .p0<-richness %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
scale_color_manual(name="treatment",
breaks=c("Acclimation", "Antibiotics", "FMT1"),
values=c("#003a45",'#008080', "#76b183")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)
p1<-neutral %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
scale_color_manual(name="treatment",
breaks=c("Acclimation", "Antibiotics", "FMT1"),
values=c("#003a45",'#008080', "#76b183")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)
p2<-phylo %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
scale_color_manual(name="treatment",
breaks=c("Acclimation", "Antibiotics", "FMT1"),
values=c("#003a45",'#008080', "#76b183")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)
4.10 Is the gut microbiota similar to the donor after FMT?
sample_metadata <- sample_metadata %>%
mutate(treatment=case_when(
treatment %in% c("FMT1", "FMT2") ~ "FMT",
TRUE ~ treatment
))
samples_to_keep <- sample_metadata %>%
filter(type == "Treatment" & treatment %in% c("FMT", "Transplant")) %>%
select(sample) %>%
pull()
subset_meta <- sample_metadata %>%
filter(type == "Treatment" & treatment %in% c("FMT", "Transplant"))4.10.1 Alpha diversity
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
mutate(metric=factor(metric,levels=c("neutral","richness","phylogenetic"))) %>%
filter(type=="Treatment" & treatment %in% c( "Transplant", "FMT")) %>%
ggplot(aes(y = value, x = treatment, color=treatment, fill=treatment)) +
geom_jitter(width = 0.2, show.legend = FALSE) +
geom_boxplot(width = 0.5, alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
scale_color_manual(name="treatment",
breaks=c("FMT","Transplant"),
values=c('#40714b', "#d5b52c")) +
scale_fill_manual(name="treatment",
breaks=c("FMT","Transplant"),
values=c('#40714b50', "#d5b52c50")) +
facet_wrap(. ~ metric, scales = "free") +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")
Richness
Modelq0GLMMNB <- glmer.nb(richness ~ treatment+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)
emmeans::emmeans(Modelq0GLMMNB, pairwise ~ treatment)Neutral
Model_neutral <- lme(fixed = neutral ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_neutral)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
173.0404 178.263 -81.5202
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 4.263247 9.322271
Fixed effects: neutral ~ treatment
Value Std.Error DF t-value p-value
(Intercept) 17.310151 3.416951 13 5.065963 0.0002
treatmentAntibiotics -8.868175 4.744329 13 -1.869216 0.0843
treatmentFMT1 7.289358 4.552796 13 1.601073 0.1334
Correlation:
(Intr) trtmnA
treatmentAntibiotics -0.596
treatmentFMT1 -0.621 0.457
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.93185401 -0.56384573 -0.04015247 0.47657219 1.55593253
Number of Observations: 24
Number of Groups: 9 Phylogenetic
Model_phylo <- lme(fixed = phylogenetic ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_phylo)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
92.7946 98.01721 -41.3973
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 0.5507824 1.405346
Fixed effects: phylogenetic ~ treatment
Value Std.Error DF t-value p-value
(Intercept) 5.488256 0.5031412 13 10.907984 0.0000
treatmentAntibiotics -1.198294 0.7137109 13 -1.678962 0.1170
treatmentFMT1 -1.351708 0.6855445 13 -1.971729 0.0703
Correlation:
(Intr) trtmnA
treatmentAntibiotics -0.611
treatmentFMT1 -0.636 0.456
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.00120146 -0.44426020 0.00785757 0.38719755 1.97397308
Number of Observations: 24
Number of Groups: 9 4.10.2 Beta diversity
Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.11100 0.111002 11.752 999 0.001 ***
Residuals 28 0.26447 0.009445
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(richness ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(richness))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(richness))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 1.184578 | 0.1548451 | 5.13002 | 0.001 |
| Residual | 28 | 6.465508 | 0.8451549 | NA | NA |
| Total | 29 | 7.650086 | 1.0000000 | NA | NA |
Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.04408 0.044082 3.1744 999 0.093 .
Residuals 28 0.38883 0.013887
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(neutral ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(neutral))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(neutral))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 1.147077 | 0.1608783 | 5.368223 | 0.001 |
| Residual | 28 | 5.983012 | 0.8391217 | NA | NA |
| Total | 29 | 7.130089 | 1.0000000 | NA | NA |
Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.00007 0.0000689 0.0044 999 0.941
Residuals 28 0.43637 0.0155847 adonis2(phylo ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(phylo))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(phylo))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.1298187 | 0.1018776 | 3.17615 | 0.001 |
| Residual | 28 | 1.1444432 | 0.8981224 | NA | NA |
| Total | 29 | 1.2742619 | 1.0000000 | NA | NA |
p0<-richness %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
scale_color_manual(name="treatment",
breaks=c("FMT","Transplant"),
values=c('#40714b', "#d5b52c")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)
p1<-neutral %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
scale_color_manual(name="treatment",
breaks=c("FMT","Transplant"),
values=c('#40714b', "#d5b52c"))+
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)
p2<-phylo %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
scale_color_manual(name="treatment",
breaks=c("FMT","Transplant"),
values=c('#40714b', "#d5b52c")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)
4.10.3 Structural zeros
FMT_samples <- sample_metadata %>%
filter(type == "Treatment" & treatment == "FMT") %>%
dplyr::select(sample) %>%
pull()
Transplant_samples <- sample_metadata %>%
filter(type == "Treatment" & treatment =="Transplant") %>%
dplyr::select(sample) %>% pull()
structural_zeros <- genome_counts_filt %>%
select(one_of(c("genome",subset_meta$sample))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
rowwise() %>% #compute for each row (genome)
mutate(all_zeros_FMT = all(c_across(all_of(FMT_samples)) == 0)) %>%
mutate(all_zeros_Transplant = all(c_across(all_of(Transplant_samples)) == 0)) %>%
mutate(average_FMT = mean(c_across(all_of(FMT_samples)), na.rm = TRUE)) %>%
mutate(average_Transplant = mean(c_across(all_of(Transplant_samples)), na.rm = TRUE)) %>%
filter(all_zeros_FMT == TRUE || all_zeros_Transplant==TRUE) %>%
mutate(present = case_when(
all_zeros_FMT & !all_zeros_Transplant ~ "Transplant",
!all_zeros_FMT & all_zeros_Transplant ~ "FMT",
!all_zeros_FMT & !all_zeros_Transplant ~ "None",
TRUE ~ NA_character_
)) %>%
mutate(average = ifelse(present == "FMT", average_FMT, average_Transplant)) %>%
dplyr::select(genome, present, average) %>%
left_join(genome_metadata, by=join_by(genome==genome)) %>%
arrange(present,-average)Structural zeros in each group
fmt <- structural_zeros %>%
filter(present=="FMT") %>%
count(phylum, name = "FMT") %>%
arrange(desc(FMT))
fmt_f <- structural_zeros %>%
filter(present=="FMT") %>%
count(family, name = "FMT") %>%
arrange(desc(FMT)) structural_zeros %>%
filter(present=="Transplant") %>%
count(phylum, name = "Transplant") %>%
arrange(desc(Transplant)) %>%
full_join(.,fmt, by="phylum" ) %>%
mutate(across(everything(), ~ replace_na(., 0))) %>%
as.data.frame() %>%
summarise(across(where(is.numeric), ~ sum(.x, na.rm = TRUE))) Transplant FMT
1 52 55Phyla to which the structural zeros belong in each group
structural_zeros %>%
filter(present=="Transplant") %>%
count(phylum, name = "Transplant") %>%
arrange(desc(Transplant)) %>%
full_join(.,fmt, by="phylum" ) %>%
mutate(across(everything(), ~ replace_na(., 0))) %>%
tt()| phylum | Transplant | FMT |
|---|---|---|
| p__Bacillota_A | 20 | 21 |
| p__Bacillota | 16 | 10 |
| p__Pseudomonadota | 6 | 11 |
| p__Bacteroidota | 3 | 6 |
| p__Desulfobacterota | 2 | 2 |
| p__Bacillota_B | 1 | 0 |
| p__Chlamydiota | 1 | 0 |
| p__Cyanobacteriota | 1 | 0 |
| p__Fusobacteriota | 1 | 0 |
| p__Verrucomicrobiota | 1 | 2 |
| p__Actinomycetota | 0 | 1 |
| p__Bacillota_C | 0 | 1 |
| p__Campylobacterota | 0 | 1 |
Families to which the structural zeros belong in each group
structural_zeros %>%
filter(present=="Transplant") %>%
count(family, name = "Transplant") %>%
arrange(desc(Transplant)) %>%
full_join(.,fmt_f, by="family" ) %>%
mutate(across(everything(), ~ replace_na(., 0))) %>%
tt()| family | Transplant | FMT |
|---|---|---|
| f__Lachnospiraceae | 7 | 9 |
| f__Erysipelotrichaceae | 6 | 5 |
| f__UBA660 | 6 | 0 |
| f__Enterobacteriaceae | 5 | 2 |
| f__CAG-508 | 3 | 0 |
| f__Ruminococcaceae | 3 | 3 |
| f__Anaerovoracaceae | 2 | 0 |
| f__Coprobacillaceae | 2 | 0 |
| f__Desulfovibrionaceae | 2 | 2 |
| f__Oscillospiraceae | 2 | 1 |
| f__Tannerellaceae | 2 | 1 |
| f__UBA1242 | 2 | 0 |
| f__ | 1 | 3 |
| f__Akkermansiaceae | 1 | 0 |
| f__CAG-239 | 1 | 2 |
| f__Chlamydiaceae | 1 | 0 |
| f__Enterococcaceae | 1 | 3 |
| f__Fusobacteriaceae | 1 | 0 |
| f__Gastranaerophilaceae | 1 | 0 |
| f__Marinifilaceae | 1 | 0 |
| f__Mycoplasmoidaceae | 1 | 1 |
| f__Peptococcaceae | 1 | 0 |
| f__Anaerotignaceae | 0 | 2 |
| f__Bacteroidaceae | 0 | 2 |
| f__LL51 | 0 | 2 |
| f__UBA3700 | 0 | 2 |
| f__Acutalibacteraceae | 0 | 1 |
| f__Arcobacteraceae | 0 | 1 |
| f__CAG-274 | 0 | 1 |
| f__CALVMC01 | 0 | 1 |
| f__Devosiaceae | 0 | 1 |
| f__Mycobacteriaceae | 0 | 1 |
| f__Pumilibacteraceae | 0 | 1 |
| f__RUG11792 | 0 | 1 |
| f__Rhizobiaceae | 0 | 1 |
| f__Rikenellaceae | 0 | 1 |
| f__Sphingobacteriaceae | 0 | 1 |
| f__Streptococcaceae | 0 | 1 |
| f__UBA1997 | 0 | 1 |
| f__UBA3830 | 0 | 1 |
| f__Weeksellaceae | 0 | 1 |
4.10.3.1 Functional capacities of the structural zeros
#Aggregate bundle-level GIFTs into the compound level
GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% structural_zeros$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>%
select_if(~ !is.numeric(.) || sum(.) != 0)
#Aggregate element-level GIFTs into the function level
GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)
#Aggregate function-level GIFTs into overall Biosynthesis, Degradation and Structural GIFTs
GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)
sample_sub <- sample_metadata %>%
mutate(treatment=case_when(
treatment %in% c("Post-FMT1", "Post-FMT2") ~ "FMT",
TRUE ~ treatment
))%>%
filter(type == "Treatment" & treatment %in% c("FMT", "Transplant"))
genome_counts_row <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>%
filter(genome %in% structural_zeros$genome) %>%
select(one_of(c("genome",sample_sub$sample))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
column_to_rownames(., "genome")
GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_row,GIFT_db)element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
significant_elements <- element_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))
element_gift_t <- element_gift %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "trait")
element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "sample")%>%
left_join(., sample_metadata[c(12,13)], by = join_by(sample == sample))
difference_table <- element_gift_filt %>%
select(-sample) %>%
group_by(treatment) %>%
summarise(across(everything(), mean)) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>%
arrange(Function) %>%
mutate(Difference=FMT-Transplant)%>%
mutate(group_color = ifelse(Difference <0, "Transplant","FMT")) difference_table %>%
ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) +
geom_col() +
# geom_point(size=4) +
#scale_fill_manual() +
geom_hline(yintercept=0) +
coord_flip()+
theme(axis.text = element_text(size = 10),
axis.title = element_text(size = 12),
legend.position = "right",
panel.background = element_blank(),
panel.grid.major = element_line(size = 0.15, linetype = 'solid',
colour = "grey"))+
xlab("Microbial Functional Capacity") +
ylab("Mean difference")+
labs(fill="Elevation")
4.10.4 Differential abundance analysis: composition
phylo_samples <- sample_metadata %>%
mutate(treatment=case_when(
treatment %in% c("FMT1", "FMT2") ~ "FMT",
TRUE ~ treatment
))%>%
filter(type == "Treatment" & treatment %in% c("FMT", "Transplant")) %>%
column_to_rownames("sample") %>%
sample_data()
phylo_genome <- genome_counts_filt %>%
filter(!genome %in% structural_zeros$genome) %>%
select(one_of(c("genome",rownames(phylo_samples)))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
column_to_rownames("genome") %>%
otu_table(., taxa_are_rows = TRUE)
phylo_taxonomy <- genome_metadata %>%
filter(genome %in% rownames(phylo_genome)) %>%
column_to_rownames("genome") %>%
dplyr::select(domain,phylum,class,order,family,genus,species) %>%
as.matrix() %>%
tax_table()
physeq_genome_filtered <- phyloseq(phylo_genome, phylo_taxonomy, phylo_samples)ancom_rand_output = ancombc2(data = physeq_genome_filtered,
assay_name = "counts",
tax_level = NULL,
fix_formula = "treatment",
p_adj_method = "holm",
pseudo_sens = TRUE,
prv_cut =0,
lib_cut = 0,
s0_perc = 0.05,
group = NULL,
struc_zero = FALSE,
neg_lb = FALSE,
alpha = 0.05,
n_cl = 2,
verbose = TRUE,
global = FALSE,
pairwise = FALSE,
dunnet = FALSE,
trend = FALSE,
iter_control = list(tol = 1e-5, max_iter = 20, verbose = FALSE),
em_control = list(tol = 1e-5, max_iter = 100),
lme_control = lme4::lmerControl(),
mdfdr_control = list(fwer_ctrl_method = "holm", B = 100),
trend_control = NULL)taxonomy <- data.frame(physeq_genome_filtered@tax_table) %>%
rownames_to_column(., "taxon") %>%
mutate_at(vars(order, phylum, family, genus, species), ~ str_replace(., "[dpcofgs]__", ""))
ancombc_rand_table_mag <- ancom_rand_output$res %>%
dplyr::select(taxon, lfc_treatmentTransplant, p_treatmentTransplant) %>%
filter(p_treatmentTransplant < 0.05) %>%
dplyr::arrange(p_treatmentTransplant) %>%
merge(., taxonomy, by="taxon") %>%
mutate_at(vars(phylum, species), ~ str_replace(., "[dpcofgs]__", ""))%>%
dplyr::arrange(lfc_treatmentTransplant)
colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
mutate_at(vars(phylum), ~ str_replace(., "[dpcofgs]__", "")) %>%
right_join(taxonomy, by=join_by(phylum == phylum)) %>%
dplyr::select(phylum, colors) %>%
mutate(colors = str_c(colors, "80")) %>% #add 80% alpha
unique() %>%
dplyr::arrange(phylum)
tax_table <- as.data.frame(unique(ancombc_rand_table_mag$phylum))
colnames(tax_table)[1] <- "phylum"
tax_color <- merge(tax_table, colors_alphabetic, by="phylum")%>%
dplyr::arrange(phylum) %>%
dplyr::select(colors) %>%
pull()Differential abundance MAGs in each group
fmt_count <- ancombc_rand_table_mag %>%
filter(lfc_treatmentTransplant<0) %>%
count(phylum, name = "FMT") %>%
arrange(desc(FMT))
ancombc_rand_table_mag %>%
filter(lfc_treatmentTransplant>0) %>%
count(phylum, name = "Transplant") %>%
arrange(desc(Transplant)) %>%
full_join(.,fmt_count, by="phylum")%>%
mutate(across(where(is.numeric), ~ replace_na(., 0))) phylum Transplant FMT
1 Bacteroidota 13 5
2 Bacillota_A 4 17
3 Bacillota 3 1
4 Campylobacterota 1 1
5 Desulfobacterota 1 2
6 Spirochaetota 1 0
7 Pseudomonadota 0 5
8 Cyanobacteriota 0 2
9 Bacillota_B 0 1
10 Bacillota_C 0 1ancombc_rand_table_mag %>%
filter(lfc_treatmentTransplant>0) %>%
count(phylum, name = "Transplant") %>%
arrange(desc(Transplant)) %>%
full_join(.,fmt_count, by="phylum")%>%
mutate(across(where(is.numeric), ~ replace_na(., 0))) %>%
as.data.frame() %>%
summarise(across(where(is.numeric), ~ sum(.x, na.rm = TRUE))) Transplant FMT
1 23 35ancombc_rand_table_mag%>%
mutate(genome=factor(taxon,levels=ancombc_rand_table_mag$taxon)) %>%
ggplot(., aes(x=lfc_treatmentTransplant, y=forcats::fct_reorder(genome,lfc_treatmentTransplant), fill=phylum)) + #forcats::fct_rev()
geom_col() +
scale_fill_manual(values=tax_color) +
geom_hline(yintercept=0) +
# coord_flip()+
theme(panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 14, face = "bold"),
legend.text = element_text(size = 12),
legend.title = element_text(size = 14, face = "bold"),
legend.position = "right", legend.box = "vertical")+
xlab("log2FoldChange") +
ylab("Species")+
guides(fill=guide_legend(title="Phylum"))
4.10.5 Differences in the functional capacities
GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts_filt$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>%
select_if(~ !is.numeric(.) || sum(.) != 0)
GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)
GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)
sample_sub <- sample_metadata %>%
mutate(treatment=case_when(
treatment %in% c("Post-FMT1", "Post-FMT2") ~ "FMT",
TRUE ~ treatment
))%>%
filter(type == "Treatment" & treatment %in% c("FMT", "Transplant"))
genome_counts_row <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>%
select(one_of(c("genome",sample_sub$sample))) %>%
column_to_rownames(., "genome")
GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_row,GIFT_db)4.10.5.1 Community elements differences:
MCI
GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
treatment MCI sd
<chr> <dbl> <dbl>
1 FMT 0.359 0.0259
2 Transplant 0.356 0.0432MCI <- GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.94871, p-value = 0.1561
Wilcoxon rank sum exact test
data: value by treatment
W = 128, p-value = 0.4826
alternative hypothesis: true location shift is not equal to 0element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
significant_elements <- element_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))
element_gift_t <- element_gift %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "trait")
element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "sample")%>%
left_join(., sample_metadata[c(12,13)], by = join_by(sample == sample))
difference_table <- element_gift_filt %>%
select(-sample) %>%
group_by(treatment) %>%
summarise(across(everything(), mean)) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>%
arrange(Function) %>%
mutate(Difference=FMT-Transplant)%>%
mutate(group_color = ifelse(Difference <0, "Transplant","FMT")) means_gift <- element_gift_filt %>%
select(-treatment) %>%
pivot_longer(!sample, names_to = "Elements", values_to = "abundance") %>%
left_join(sample_metadata, by=join_by(sample==sample)) %>%
group_by(treatment, Elements) %>%
summarise(mean=mean(abundance))
log_fold <- means_gift %>%
group_by(Elements) %>%
summarise(
logfc_fmt_transplant = log2(mean[treatment == "FMT"] / mean[treatment == "Transplant"])
) %>%
left_join(., difference_table, by="Elements")difference_table %>%
ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) +
geom_col() +
# geom_point(size=4) +
scale_fill_manual(values=c('#40714b', "#d5b52c")) +
geom_hline(yintercept=0) +
coord_flip()+
theme(axis.text = element_text(size = 10),
axis.title = element_text(size = 12),
legend.position = "right",
panel.background = element_blank(),
panel.grid.major = element_line(size = 0.15, linetype = 'solid',
colour = "grey"))+
xlab("Microbial Functional Capacity") +
ylab("Mean difference")+
labs(fill="Treatment")
log_fold %>%
filter(Elements!="D0611") %>%
ggplot(aes(x=forcats::fct_reorder(Function,logfc_fmt_transplant), y=logfc_fmt_transplant, fill=group_color)) +
geom_col() +
#scale_fill_manual() +
geom_hline(yintercept=0) +
coord_flip()+
theme(axis.text = element_text(size = 10),
axis.title = element_text(size = 12),
legend.position = "right",
panel.background = element_blank(),
panel.grid.major = element_line(size = 0.15, linetype = 'solid',
colour = "grey"))+
xlab("Microbial Functional Capacity") +
ylab("Log_fold")+
labs(fill="Treatment")
#### Community functions differences
MCI
GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
treatment MCI sd
<chr> <dbl> <dbl>
1 FMT 0.349 0.0221
2 Transplant 0.355 0.0377MCI <- GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.87044, p-value = 0.001713
Wilcoxon rank sum exact test
data: value by treatment
W = 121, p-value = 0.6801
alternative hypothesis: true location shift is not equal to 0function_gift <- GIFTs_functions_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
merge(., sample_metadata[c(12,13)], by="sample")unique_funct_db<- GIFT_db[c(3,4,5)] %>%
distinct(Code_function, .keep_all = TRUE)
significant_functional <- function_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
left_join(.,unique_funct_db[c(1,3)],by = join_by(trait == Code_function))
significant_functional# A tibble: 3 × 4
trait p_value p_adjust Function
<chr> <dbl> <dbl> <chr>
1 B04 0.000220 0.00441 SCFA biosynthesis
2 B10 0.00284 0.0189 Antibiotic biosynthesis
3 S02 0.00174 0.0174 Appendages